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AFIT/GEP/PH/80D-7 


Abstract 


The  mode  eigenvalue  equation  for  an  unstable  strip 
laser  resonator  is  developed  from  scalar  diffraction 
theory.  The  field  distributions  are  expressed  as  a  series 
and  the  integral  is  then  evaluated  using  a  first  order 
approximation  to  the  method  of  stationary  phase.  The  re¬ 
sulting  approximate  closed  form  is  rearranged  to  form  an 
eigenvalue  polynomial,  the  roots  of  which  are  the  mode 
eigenvalues.  Eigenfunction  expressions  are  then  developed 
using  a  second  order  approximation  to  the  method  of  station¬ 
ary  phase.  Modifications  to  these  expressions  are  then 
made  to  account  for  the  presence  of  uniform  gain  in  the  res¬ 
onator. 

The  results  of  a  computer  program  using  the  derived 
expressions  are  presented.  Comparisons  to  previously  pub¬ 
lished  results  are  made  for  the  bare  cavity  case,  and  re¬ 
sults  for  the  loaded  cavity  case  follow. 
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ANALYSIS  OF  MOOES 


IN  AN  UNSTABI  E 
STRIP  LASER  RESONATOR 

I.  Introduction 


Background 

An  unstable  laser  resonator  is  a  resonator  in  which  the 
geometric  path  of  a  paraxial  ray  traveling  back  and  forth  be¬ 
tween  the  two  mirrors  is  unbounded  in  an  infinite  number  of 
passes.  This  is  opposed  to  a  stable  resonator,  in  which  the 
ray  path  is  bounded.  Any  ray  inside  an  unstable  resonator 
will  eventually  take  on  a  direction  from  which  it  will  not 
come  into  contact  with  either  mirror,  and  thus  leave  the 
cavity.  In  this  type  of  resonator,  the  product  of  the  res¬ 
onator  mirror  g  parameters,  where 


lies  outside  the  stable  range  of 

0  <  9^2  *  1 

The  utility  or  benefits  of  unstable  resonators,  for  in¬ 
stance  large  mode  volume  and  minimally  transmitting  optics 
(Ref  10:353),  require  that  some  method  of  mode  analysis  be 
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available.  Several  methods  are  available,  but  have  various 
drawbacks,  such  as  excessive  computer  processor  time  re¬ 
quirements,  or  limited  applicability. 

Horwitz  (Ref  6)  developed  a  method  whereby  the  mode 
eigenvalue  equation  for  an  unstable  strip  resonator,  mod¬ 
ified  from  the  original,  developed  by  Fox  and  Li  (Ref  4), 
was  simplified  by  using  first  a  series  of  functions  found 
through  asymptoti c  analysis  to  approximate  the  field  in 
the  resonator,  and  then  the  method  of  stationary  phase  to 
approximate  the  integral.  Butts  and  Avizonis  (Ref  2)  clar¬ 
ified  this  approach  and  modified  it  to  allov.'  consideration 
of  a  resonator  with  circular  mirrors.  However,  neither 
allowed  for  the  inclusion  of  a  gain  medium  in  the  cavity. 

Ob jecti ves 

The  objective  of  this  thesis  is  to  develop  a  computer 
code  allowing  analysis  of  modes  in  an  unstable  resonator 
and  to  then  utilize  that  code  in  performing  said  analysis. 

The  code  is  to  be  developed  for  a  strip  resonator  and  ac¬ 
count  for  both  bare  and  loaded  cavity  cases. 

Assumptions 

To  facilitate  modeling  of  the  unstable  resonator,  cer¬ 
tain  simplifying  assumptions  will  be  made: 

1.  Scalar  diffraction  theory  will  be  used  to  describe 
the  physical  situation  in  the  resonator.  This  is  reasonable. 
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since  the  dimensions  of  laser  resonators  are  large  compared 
to  optical  wavelengths. 


2 .  The  Fresnel  approximation  to  the  Kirchoff-Fresnel 
formula  is  valid.  Resonator  cavity  lengths  make  this  an 
acceptable  assumption. 

3.  In  a  Cartesian  system,  diffraction  integrals  and 
mode  eigenfunctions  are  separable.  This  allows  a  1-D  strip 
resonator  to  be  utilized  in  the  following  development. 

4.  One  of  the  resonator  mirrors  is  very  much  larger 
than  the  beam  width  on  that  mirror.  In  other  words,  that  the 
height  of  this  mirror  be  considered  infinite.  This  is  not 

an  impossible  physical  constraint. 

5.  The  modes  in  the  strip  resonator  consist  of  a  fun¬ 
damental  cylindrical  wave  modified  by  a  finite  number  of 
diffraction  effects.  This  assumption  is  supported  by  early 
analysis  of  unstable  resonators.  (Ref  9:279) 


Procedure 

This  thesis  will  start  with  the  Kirchoff-Fresnel  dif¬ 
fraction  formula  and  develop,  following  Horwitz  (Ref  6:1529), 
the  eigenvalue  equation  for  a  strip  resonator 


Ag(x) 


(1.3.1) 
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where  the  eigenfunctions  g(x)  are  the  weighting  functions 
of  the  basic  cylindrical  wave  assumed  to  be  present  in  the 
resonator,  that  is,  if  the  total  field  is  described  by 
u(x)  then  (Ref  11 ) , 


-  iirN^x2 

u(x)  =  g(x)e 


(1.3.2) 


-inN^x2 

e  being  the  phase  curvature  term  of  the  basis  cylin 

drical  wave. 

The  eigenvalues  will  be  found  by  developing  a  suitable 
relation  from  the  eigenvalue  equation.  The  total  field  in 
the  cavity  is  first  assumed  to  consist  of  a  unit  amplitude 
cylindrical  wave  plus  a  finite  series  of  diffraction  sup¬ 
plements.  This  is  stated  in  terms  of  g(x)  as  (Ref  2) 


g(x) 


+ 


N 


X  cnHn(x) 
n  =  1 


(1.3.3) 


This  expression  is  substituted  in  the  eigenvalue 
equation  and  then  an  approximation  to  the  integral  is  devel 
oped  using  a  first  approximation  to  the  method  of  station¬ 
ary  phase.  The  resulting  relation  will  allow  the  eigen¬ 
values  to  be  expressed  as  roots  of  a  polynomial  with  deter¬ 
minable  coefficients.  The  roots  can  be  found  by  using  a 
general  root-finding  routine. 

An  eigenfunction  expression  may  then  be  found  by  using 
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the  original  assumption  ie.,  equation  (1.3.3).  However, 
inherent  limitations  of  the  first  stationary  phase  approx¬ 
imation  confining  applicable  x  values  in  this  relation 
require  the  development  of  a  better  approximation  to  the 
integral.  The  higher  order  expression  will  be  developed, 
using  the  higher  order  approximation  to  the  method  of 
stationary  phase  (Ref  1),  enabling  the  evaluation  of  the 
eigenfunctions  throughout  a  continuous  range  of  x  values. 

This  thesis  will  then  seek  to  modify  the  bare  cavity 
expressions  to  account  for  a  gain  medium  in  the  resonator 
by  introducing  a  gain  factor,  e^*-  t  into  the  integral 
abd  by  relaxing  the  unit  amplitude  requirement  on  the  fun¬ 
damental  cylindrical  wave. 


Organi zati on 

The  derivation  of  the  basic  resonator  eigenvalue 
equation  will  be  covered  in  Chapter  II.  Chapter  III  will 
present  the  two  applications  of  that  equation:  calculation 
of  eigenvalues  and  evaluation  of  eigenfunctions.  Inclusion 
of  gain  considerations  will  be  covered  in  Chapter  IV  and 
Chapter  V  will  contain  results  of  the  computer  code.  Chap¬ 
ter  VI  will  include  conclusions  and  further  recommendations 


1 1 .  Development  of  the  Eigenvalue  Equation 


Chapter  II  addresses  the  problem  of  applying  the  Kir- 
chof f-Fresnel  diffraction  formula  to  the  desired  case  of  an 
unstable  optical  resonator.  The  development  follows  that 
in  Reference  6. 

A  steady  state  mode  will  exist  in  a  resonator  when  the 
field  value  on  one  mirror  resulting  from  one  round  trip 
through  the  resonator  multiplied  by  some  complex  constant 
is  equal  to  the  original  field  value  on  that  mirror.  Math¬ 
ematically  this  can  be  stated  as 

YE'(x.y)  =  E(x,y)  2.1.1 


where  E  is  the  original  field  distribution  on  M2  ,  the 
second  mirror,  E'  is  the  distribution  after  one  round  trip, 
and  y  is  the  constant,  in  general  complex. 

Wave  propagation  through  the  resonator  can  be  express¬ 
ed  using  scalar  diffraction  theory.  Wave  propagation  from 
a  rectangular  aperture,  dimensions  2a  x  2c  ,  on  one  plane 
to  another  plane  a  distance  L  away ,  as  seen  in  Fi g .  1 ,  is 
given  in  the  Fresnel  approximation  by 


E ( x 2 ,y2  ) 


AL 


c  a 

J  J  E ( x  i  ,y  i  )e 

-c  -a 


4t-[(xi-x2  )2  +  (yi-y2  )2] 

dx i dy i 
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2.1.2 


Figure  1 


■  Since  this  is  presented  in  a  Cartesian  system,  the 
field  distribution  can  be  separated.  This  is  done  by  as 
suming  (Ref  4:485-486), 

E(x,y)  =  U (x)U(y )  2.1 

Substitution  of  2.1.3  into  the  diffraction  formula 
yields  two  independent  diffraction  formulae. 


U(y2 ) 


i  kL  c 

e  2  J  U  (y  x ) 


"IPyi~y2 ) 


dy  i 


-c 


2.1.5 


Consideration  of  only  one  of  these  formulae  is  equiva¬ 
lent  to  considering  diffraction  from  a  strip  aperture.  No 
generality  is  lost,  however,  since  the  effects  of  a  finite 
aperture  can  be  found  from  the  product  of  two  separate  strip 
cases.  Thus  the  one  remaining  equation  is 


U  (  X  2  ) 


/ 


~2T(xi-x2 ) 


2 

U ( x 1 )dxi 


2.1.6 


Equation  2.1.6  represents  propagation  from  one  plane  to 
another.  For  this  to  correctly  represent  propagation  in  a 
resonator,  the  phase  lag  introduced  by  mirror  curvature  must 
be  accounted  for. 

The  phase  lag  introduced  by  the  mirrors  can  be  express¬ 
ed  as  a  function  of  distance  from  the  optic  axis.  This  ex¬ 
pression  can  be  derived  from  the  paraxial  lens  thickness  func¬ 
tion,  (Ref  5:80),  which  is 


A(x,y)  =  A0  -  (]T7  "  2.i .7 

Where  A  =  thickness,  A0=  maximum  lens  thickness,  x  and  y 
are  coordinates  of  the  point  where  the  ray  of  interest  is 
incident  on  the  lens,  and  Rj  and  R2  are  the  radii  of 


8 


Figure  2 


curvature  of  the  lens'  surface.  In  the  case  of  a  strip 
mirror  the  lens  equation  applys  if 

R2  =  CO 

y  =  0 

A  o  =  0 


giving 


A  (  x  )  = 


2R. 


2.1. 


for  the  ith  mirror.  The  phase  lag  at  some  particular  dis 
tance  .from  the  optic  axis  x  is  given  by 
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x  rr  -  ^  r  -  Tr[9  ixi2+92X22-2x1x2] 

U(x2)  =  e  2  J  U(Xl)e  2L  dxx 


2.1.13 


This  expression,  now  modified  to  describe  propagation 
of  U(xi)  from  Mx  to  M2  ,  can  be  used  to  set  up  two 
equations:  one  for  propagation  from  M2  to  Mx  and  the 


other  for  propagation  from  Mx  to  M2 


Combination  of 


the  two  will  then  yield  an  expression  describing  propaga¬ 
tion  of  a  field  through  one  round  trip  in  the  resonator. 
The  one  way  formulae  are 


U(x2)  =  e 


V  -fr[gixr2+g2x22-2x:x2] 

f  U  (  x  x '  )e  2L  dxr 


2.1.14 


U(Xl)  = 


rr  V  4r[g.*i2+g2*i 

-  a  2 


- 2x  !  x 2 ] 


If  2.1.15  is  substituted  for  u(xj)  in  2.1.14,  the  re¬ 
sultant  expression  will  give  the  field  on  M2  due  to  the 
propagation  of  an  original  field  on  M2  through  one  round 
trip  in  the  resonator.  Substitution  gives 


U(x2) 


ai 

e"lkL  /T/TT  j  /T/lL 
-a2 


a  i 

j  U(x£) 

-a  i 


-  4f[gixr2+g2x£2-2xrx2]  -  if[ gixr2+g2xl-2xj'x2] 

e  dx£e  ^  dxT 


2.1.16 

In  the  case  considered  in  Fig.  2,  the  assumpti on  that 
M!  is  much  bigger  than  the  beam  width  on  that  mirror  for 
any  laser  mode  that  is  likely  to  resonate,  allows  ai  to 
be  thought  of  as  essentially  infinite.  Then  2.1.16  be¬ 
comes 


U(x2) 


eikL 


a2  ® 

f  fi r 


-  a  2  -® 


|£[gixr2  +  g2x22-2xfx2'j 


-^■[g1xj2  +  g2x22-2^x2] 

e  U(x2)dx^dx2  2.1.17 


This  expression  can  be  simplified  by  extracting  the  in¬ 
terior  integral 

i  7  ■|f^9lXj’2  +  g2X22"2x^x^ 

XL  J  e 

-00 

-yr[g iXi2+g2x22-2xix2] 

e  dxj  2.1.18 
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Evaluation  of  this  integral  in  Appendix  B  yields 


rr~  -2ir?rt(29ig2-1)(x£2+X22)-2X2X^  9  .  . 

"y  2LXg  6 

Substitution  of  this  for  the  complete  kernel  in  2.1.17  in 
turn  yields 


U(x2) 


-Yryr-[(2g1g2-1) (X22+X22 )-2x2X2] 
e  Xu2gi  U(xJ)  d x 2 


2.1.20 


To  simplify  this  further,  the  definitions 


2  g  i  g  2  - 1  =  g 


2.1.21 


and 


a  2 2  _  F  2  = 

2g i XL  2  g  i 


2.1.22 


are  introduced. 

Here,  gt  and  g2  are  the  familiar  g  parameters 
and  F2  is  the  ordinary  Fresnel  number  of  the  smaller  feed¬ 
back  mirror.  The  ordinary  Fresnel  number  is  defined  as  the 
additional  length  per  pass  in  half  wavelengths  for  a  ray 
traveling  from  one  mirror's  center  to  the  other  mirror's 
edge,  compared  to  one  traveling  from  mirror  center  to  mirror 
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center.  (Ref  11:159-161). 

The  dimensions  of  the  quantities  are  also  scaled  such 
that  a2  =  l  •  These  modifications  yield 


1 

U(x2)  =  e“ikL/TF  J U  (  x  2 ) e 

-1 


-i7Tp[g(x22  +  x2  )  -2x2x2] 


dxi 


1.23 


Imposition  of  the  reproducobi 1 i ty  constraint,  equation 

~  i  k  L 

2.1.1,  and  absorption  of  the  constant  e  into  y 

yields 

1 

yU(  xz  )  =  /TF  f  U(xne-(,FC9(x^tx"2)-2^xhdx|  2.1.24 

-1 


Introducing  the  dummy  variable  y  and  dropping  sub¬ 
scripts  and  superscripts  yields 

1 

yU  (  x  )  -  /TF  f  U(y)e-i"Ft9(x2+y2)-2x>']dy  2.1.25 

-1 


To  further  simplify  this  equation,  the  following 
quantities  are  defined 


Nf  =  £(m  -  i) 


2.1.26 


and  g(x)  such  that 


U  ( x ) 


-  i  TT  N  ^  x 2 
e  g  ( x ) 


2.1.27 
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is  the  equivalent  Fresnel  number  of  the  resonator. 
The  equivalent  Fresnel  number  can  be  interpreted  as  the 
additional  path  length  per  pass  in  half  wavelengths  for  a 
ray  traveling  from  a  mirror's  virtual  center  to  the  edge 
of  the  next  mirror,  as  opposed  to  a  ray  traveling  from  the 
virtual  center  of  one  mirror  to  the  actual  center  of  the 
next  (Ref  11:159-161).  The  virtual  center  is  defined  as 
that  point  from  which  a  cylindrical  wave  would  eminate  if 
that  wave  were  to  be  reflected  from  a  feedback  mirror,  and 
then  return  to  the  original  mirror  in  the  same  form  as  when 
it  left  (Ref  9:279-280).  That  cylindrical  wave  is  then  as¬ 
sumed  to  take  the  form 


-iirN^x2 


2.1.28 


and  the  entire  wave  function  is  assumed  to  be  based  on  that 
wave,  stated  by  2.1.27.  Substitution  of  2.1.26  and  2.1.27 
into  2.1.25  yields 


-mj( 

yg(x)e 


m2  - 1 
m 


1 

/TF  J  g(y)e 


-iirFy2  ( 


m2  - 1  v 

2m  ' 


-1 


-iirF[g(x2+y2  )-2xy] 
e  dy  2.1.29 


After  some  manipulation,  detailed  in  Appendix  C, 
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2.1.29  simplifies  to 
gral  equation 


the  final  form  of  the  resonator  inte- 


where 


1 

yg ( x )  =  /i  t/irm  J  g(y)e 

-1 


1  -i  t(y-— )2 

xg(x)  =  /TUT  J  e  m  gCy)  dy 

-1 


2.1.30 


2.1.31 


2.1.32 


2.1.33 
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III.  Determination  of  Eigenvalues  and 


Evaluation  of  Eigenfunctions 


Chapter  III  is  concerned  with  solving  the  resonator 
mode  eigenvalue  equation  and  with  developing  expressions 
for  the  resulting  eigenfunctions.  The  eigenfunctions  are 
most  desireable  since  they  will  ultimately  express  field 
values  across  the  output  mirror  plane. 


where  g(x)  is  the  quantity  multiplying  the  primary  cylin¬ 
drical  wave  expressed  as  a  function  of  normal  distance  from 
the  optic  axis. 

Now  it  is  assumed  that  the  field  on  the  mirror  before 
the  round  trip,  U(y)  ,  consists  of  a  unit  amplitude  cylin¬ 
drical  wave  plus  an  infinite  series  of  edge  diffracted  waves 
given  by  some  functions  Hn(y)  (Ref"  2).  In  terms  of  g(y) 
this  is  stated  as 

oo 

g(y)  =  1  +  2  cnHn(y)  3.1.2 

n=  1 
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The  physical  basis  for  this  assumpti on  is  that  the 
original  field  on  M2  will  consist  of  that  primary  cylin¬ 
drical  wave  which  makes  the  round  trip  unchanged  plus 
other  contributions  which  are  the  diffraction  additions  to 
that  wave  from  previous  reflections.  To  make  this  viable, 
however,  it  is  then  assumed  that  the  series  terminates  when 
eventually  some  function  H^(y)  is  the  last  contribution 
that  has  any  new  effect  on  the  field,  or  that  HN+1(y)  is 
constant.  If  the  resonator  is  thought  of  as  an  infinite 
lens  train,  the  mode  components  between  the  last  two  lenses 
will  consist  of  the  basic  cylindrical  wave  and  one  diffrac¬ 
tion  effected  wave  from  each  preceeding  lens  group.  The 
series  terminates  when  the  consideration  of  another  lens 
group,  farther  back,  adds  no  more  new  information  to  the 
final  mode.  Then  the  addition  of  one  more  diffraction 
effected  wave  would  add  only  to  amplitude,  and  not  change 

the  shape  of  the  total  wave.  3.1.2  then  becomes 
N 

g(y)=l+ni1  cnHn(y)  .  A  good  approximation  is  to  let  (Ref. 
6:1533) 

N  ^  l?— !50-N-f-  3.1.3 

1  n  m 


and  the  quality  of  this  approximation  is  displayed  in 
Appendix  E. 

When  3.1.2  is  substituted  into  3.1.1,  the  result  is 


Ag(x)  =  /i  t / tt 


1  r  N 

/  -Da 

-1  l  n=  1 


(y) 


“ i t(y-m) 2 

e  m  dy 


3.1.4 
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Some  method  of  approximating  this  integral  is  needed. 
The  method  chosen  is  the  method  of  stationary  phase.  This 
method  states  that  an  integral  of  the  form 

b 

f  e'ltp(y)q(y)  dy  3.1.5 

a 

can,  when  t  is  large  and  q ( y )  is  slowly  varying,  be  ex¬ 
pressed  as  a  series,  the  first  two  terms  of  which  are  ap¬ 
proximately  (Ref .2: 1073) . 


-iir/4  i  x  -  i  t p  ( y  o )  /  2tt 

q(yo)e 


+1  q(b)  --itp(b)  q ( a )  -itp(a) 
+  t  pTFT  6  "  pT5T  6 


where  y0  is  the  point  of  stationary  phase,  ie. 

p'(yo)  =  0  .  3.1.7 


To  utilize  this  however  some  explicit  form  of  Hn(y) 
is  needed.  The  form  used  here  is  the  same  as  that  develop¬ 
ed  by  Horwitz  through  asymptotic  analysis  of  the  resonator 
integral,  2.1.33.  The  form  is  as  follows: 

Given  the  functions  (Ref  3) 


F(x,t) 


1  e  -  i  t  ( 1  -  x  )  2 

2/TTt  1  -  x 


3.1.8 
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3.1.9 


functions: 


3.1.15 


cnHn(x>  ’ 


a  F  (x)+b  G  (x) 
n  iv  n  n  ' 


Since  the  cavity  under  consideration  here  is  centered 
on  the  optic  axis,  symmetry  dictates  either  odd  or  even  field 
functions.  To  get  an  even  field  function,  then  it  is  assum¬ 
ed  that  Hn(x)  is  even.  Odd  would  require  that  Hn(x)  be 

odd . 


It  is  seen  from  3.13-3.15  that  Hn(x)  can  be  made  even 


if  a  =b 
n  n 


So,  if 


a  =  b 
n  n 


3.1.16 


then  with 


and 


However,  since 

t 


c„ 3 

an  3  bn 

3.1.17 

cnHn(x) 

'  cr,(Fn(x)+Gn(x)) 

3.1.18 

cnHn(-x) 

3  cn(Fn<-x)+Gn(-x)) 

3.1.19 

F„(x) 

3  Gn(-X> 

3.1.20 

cnHn(-x) 

"  cn(Gn(x)*Fn(x)) 

3.1.21 

3  cnHn(x) 

3.1.22 
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Thus,  Hn(x)  i s  an  even  function. 

Similarly,  Hn(x)  can  be  ^ade  odd  by  assuming 

a  =  -b  3.1.23 

n  n 

and  that 


an  =  "b  =  Cn  3.1.24 

n  n  n 

it  is  seen  that 


cnHn(x)  ’  cn(Fn(x)-Gn(x)> 


’  cn(Fnl-x)-Gn(-x> 

3.1.26 

'  cr,(Gn(x)-Fn(x>> 

3.1.27 

*  -cnH„(x> 

3.1.28 

Thus  Hn(x)  i s  an  odd  function.  One  additional  assumption 
is  that  in  the  odd  case,  the  amplitude  of  the  cylindrical 
wave  is  zero.  This  is  necessary  for  the  field  function  to 
be  odd. 

In  the  following  development,  the  even  parity  case 
will  be  the  one  dealt  with.  The  odd  parity  equations  can 
be  found  from  those  for  the  even  case  by  deleting  the  lead¬ 


ing  term  in  eq.  3.1.2  and  following  the  procedure  as  above. 
Therefore,  the  eigenvalue  equation  to  be  solved  is 


Ag(x)  =  /"i  t /tt 


l 

l 


X\2 


■itfy.i) 


1+£(anFn(y)tbnGn<y>> 


n=  1 


dy 


3.1.29 


Substitution  of  the  actual  forms  of  the  functions  allows 
explicit  forms  of  p(y)  and  g ( y )  found  in  3.1.4  and 
3.1.5  to  be  found.  Employing  equation  3.1.5,  according  to 
Appendix  A  allows  the  following  first  order  approximation 

Ag(x)  =  1+Fj ( x ) +Gj (x) 

N 

*  £(anFn+l(x)+bnGnU(x)) 
n  =  l 

N 

*F><x>I>nF„<1>+b„Gn<1>> 
n  =  1 

N 

+G1(x)^(anFn(-l)  +  bnGn(-l))  3.1.30 

n  =  l 

N  *  N 

“  1+Hi(x)+5^cnHn+1(x)+H1(x)^cnHn(l)  3.1.31 

n=l  n=l 

N 

=  *(1+]T/nHn(x))  3.1.32 

n=  1 
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24 


and  generally  it  is  seen  that 


cn+l 


3.2.5 


In  other  words 


CiX 


V 


3.2.6 


=  cNx 


3.2.7 


This  in  turn  implies  that 


c 


n 


CNX 


N-n 


3.2.8 


Equating  coefficients  of  H  (x)  now  yields 


X  C  j 


N 


1+EcnHn(1) 

n  =  l 


3.2.9 


Substituting  for  cn  and  C]  according  to  3.2.8,  gives 


N 

XcNXN_1  =  l+23cNxN’"Hn(l)  3.2.10 

n=l 


Equating  constant  terms  in  3.2.1  shows  that 

X  =  1+cnHn+1  3.2.11 
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r 


X-1  =  cnhn+i 


3.2.12 


X-l 

HN+1 


3.2.13 


Substituting  for  in 


A  (  A  -  1 )  A 
HN+1 


N-l 


3.2.10  according  to  3.2.13  yields 

1  +  Y  (A~1)xN-n-  H  (1)  3.2.14 

fa  N+l  n 


or 

N 

AN(A-1)  =  Hn+1  +  ( A- 1 )  ^XN"nHn(l)  3.2.15 

n  =  l 

which  is  a  polynomial  in  the  complex  variable  A  .  Its 
roots  can  be  determined  from  any  root-finding  subroutine, 
since  its  coefficients  all  involve  known  quantities  such  as 

H„<i) 

or  the  constant 

HN+1 

It  is  from  this  polynomial  that  the  mode  eigenvalues  of  the 
resonator  are  determined.  A  preliminary  evaluation  of  the 
eigenfunction  for  a  particular  mode  can  be  made  by  substi¬ 
tuting  into  equation  3.1.2  the  values  for  cn  ,  which  are 
given  by 
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cn  '  CNX 


N- n 


( X- 1 )  ,N-n 


3.2.16 


However,  due  to  the  singularities  in  the  first  approximation 
to  the  integral,  3.1.6,  whenever  x  approaches  y0  ,  this 
particular  expression  for  the  eigenfunction  is  invalid.  This 
problem  will  be  remedied  in  the  next  section. 

The  odd  parity  solution  is  given  by  3.1.33,  and  the  poly¬ 
nomial  development  for  that  case  is  as  follows. 

After  equating  the  coefficients  of  Hn(x),  n  t  1 ,  it  is 
seen  that  the  same  relations  arise  as  3.2.2  -  3.2.8. 

Equating  coefficients  of  Hj(x)  indicates  that 

N 

Xc,  =  5^cnHn(l)  3.2.17 

n  =  l 

Equating  constant  terms  indicates  that 


0 


CNHN+1 


This  is  only  reasonable  since  the  condition  imposed  on  N 
namely  that  is  a  constant,  also  implies  that 


FN+1  =  Constant  =  gn+i 


3.2.18 
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4 

i 


i 


and  since 


HN+1  =  FN+1  ‘  GN+1 


3.2 


HN+1  =  0  3-2 

This  indicates  that  cn  is  completely  arbitrary  since 
are  no  other  restrictions  imposed  by  either  3.2.8  or  3. 
If  cn  is  indeed  arbitrary,  and 


then  c  can  be  chosen  such  that 

n 


V  •  vN  ■  1 


leaving  the  relation 


3.2 


3.2 


which  can  be  used  in  the  limited  range  eigenfunction  ex 
pression  for  the  odd  parity  case. 

The  polynomial  is  developed  by  substitution  for  c 
of3.2.8in3.2.17giving 

N 

ic’=ECnHn(1)  3-2 

n=  1 


.19 

.20 

there 

2.17. 

.21 

.22 

.23 

n 

.24 
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3.2.25 


xcNx 


-‘■Ev'-Vi 

n*  1 


A 


N 


N 


I>N-\<n 


3.2.26 


Development  of  Eigenfunction 
Expressions  Valid  for  All  X 

To  develop  an  eigenfunction  expression  valid  for  all 
x  ,  it  is  first  necessary  to  return  to  the  original  equa¬ 
tion,  2.1.32,  which  is 


Ag(x) 


-  i 

g(y)e 


dy 


3.3.1 


Since  the  eigenvalues  are  known  or  can  be  determined, 
it  can  then  be  said  that 


g(x) 


i/Tt77 


t(y-i)’ 


dy 


3.3.2 


One  might  question  the  validity  of  this  expression,  since 
A  was  determined  from  the  first  order  approximation.  How¬ 
ever,  that  previous  approxima ti on  yields  perfectly  valid 
values  for  A  ,  because  all  that  determines  the  mode 
eigenvalue  is  the  field  on  the  smaller  feedback  mirror, 
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where  x£l  .  In  this  region,  the  approximation  is  always 
valid.  Therefore  the  A ’ s  are  perfectly  valid. 

All  expressions  and  quantities  on  the  right  side  of 
3.3.1  are  known,  and  therefore  one  can  once  again  utilize 
the  method  of  stationary  phase,  but  in  a  second  approxima¬ 
tion,  yielding  an  expression  no  longer  as  simple  as  3.1.6 
but  one  that  is  valid  for  all  x  .  The  higher  approximation 
to  the  integral  is  given  by  (Ref  1) 


y0 


vnen 


is  such  that 


when  y9  is  such  that 


when  y0  is  such  that 


Here,  E*  is  the  complex  conjugate  of  the  Fresnel  integral. 
In  both  the  even  and  the  odd  cases  the  integral 


J-  /?  t/TT 


i 

/ 


^  1  2 


£< 

n=  1 


a„Fn<*>tbnG„<*>> 


dy 


3.3.9 


must  be  evaluated.  Specific  differences  for  even  and  odd 
cases  will  be  treated  later.  Manipulating  3.3.5  yields 


(a  F  (y)+b  G  ( y ) }  dy 
vnnw/  nnw  J 


3.3.10 


To  make  use  of  equations  3.3.3,  3.3.5  and  3.3.7,  it  is 

f"  h 

necessary  to  get  p(y)  and  q(y)  expressions  for  the  nl 
term  in  the  series.  Substitution  of  the  explicit  forms  of 
the  Fn  and  Gn  functions  yields  an  integral  of  the  form 


x 

l 


n  - 1 


2 /iirt 


1  -  X- 

1  mn 


V 


bj/TT 

n*  n -  1  e 


2/TiT 


1  +  X- 
1  mn 


dy 


3.3.11 


where  has  been  previously  defined  in  3.1.12. 

Upon  consideration  of  the  term  involving  the  an  con¬ 
stants,  it  is  seen  that 
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3.3.12 


p(y)  = 


(y--)  + 

m  m 


d-£rr) 


n- 1 


P'(y)  =  2(y-£)  -  ~tt  d-^n) 


P~(y)  =  2 


Vlm 


2n 


Ttr 

1  mn 


Solving  for  y0  yields 


..  x  1  +  y  o  _  n 

'  m  '  nVlm  \.^Zn 


yod+^ — ^rr)  -  £ 


X  +  l 


mn_  ^mn '  m  mn-lmn 


3.3.13 


3.3.14 


3.3.15 


3.3.16 


3.3.17 


3.3.19 


p(y)  =  ( y -— )  + 


( 1  +  &r) 


n- 1 


P'(y)  =  2(y-J)  +  pr  (l+Jrr) 


p""(y )  =  2  +  — 2^- 


m  m 


n-  1 


q(y)  = 


l  +  X— 
1  mn 


Solving  for  y0  ,  it  is  seen  that 


These  expressions  can  now  be  substituted 
all  approximations  to  the  integral.  However, 
careful  consideration  of  the  y0  values  must 


3.3.20 

3.3.21 

3.3.22 


3.3.23 

3.3.24 

3.3.25 

into  the  over 
in  evaluation 
be  taken,  in 
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order  that  the  proper  form  of  the  approximation  is  used 

This  is  rather  complicated,  since  there  are  two  yo's 

one  for  the  a  term,  and  one  for  the  b  term. 

n  n 


In  order  to  simplify  things,  let 


=  ATERM 


3.3 


and  for  the 

_  bn Vmn-1 
2/vttT 


3.3.32 


These  expressions  are  added  to  ATERM  or  BTERM,  which¬ 
ever  is  required.  In  this  way  the  complete  expression  for 
the  integral 

N  1 

lE  /'  ta„Fn(y)+b„Gn(y))  dy  3.3.33 

n  =  l  =1 


can  be  stated. 

However,  in  the  even  parity  case,  one  more  modification 
must  be  made.  The  term  involving  the  *  1  *  must  be  added  to 
the  expression.  Explicity,  that  term  is 


1 

A 


t(y-i) 


dy 


3.3.34 


C 


From  this  it  is  seen  that 

p(y)  -  (y  -  £)! 
p'<y)  -  z(y  -  £> 


3.3.35 


3.3.36 
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a 


p"(y) 

q(y) 


and 


y<>  =  - 


3.3.37 

3.3.38 


3.3.39 


Substituting  into  3.3.2  for  y0  ^  -1  gives 

x  \  2  i  /  _ 


x  /TT7T 


e  e 


-it(-i-i) 


XX2 


-e 


i t/4*  4( -  1-4) 2 


nr 


E* 


2tt 


2(  i~) 
m 


1-i 


3.3.40 


If  y0  -  1  ,  once  again,  all  that  needs  to  be  done  is  to 
change  the  sign  of  the  terms.  If  y0  is  such  that 

3.3.30  is  satisfied,  then  only  the  second  -  term  is 
changed  in  sign  and  the  stationary  phase  point  contribution 
term  is  added.  That  term  is  given  by 


r1*/4  /77T 


3.3.41 


Thus,  the  higher  order  approximation  expressions  for 
the  eigenfunctions 


x  /TT7¥ 


N  1  -u(*-i)! 

e  m  (1+anFn(y)+bnGn(y))  dy  3.3.42 

n  =  1  -1 


N  1 

£  / 
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I 


in  the  even  case,  and 


i  /TT7?  £  /  e  (anFn(y)+bnGn(y))  dy 

n  =  l  -1 

in  the  odd  case  are  expressed,  per  3. 3. 2-3. 3. 4. 
pressions  are  valid  for  all  x  .  In  this  way, 
across  the  output  mirror  plane  can  be  evaluated 
lows  that  intensities  are  then  given  by 

I  =  £*E  =  g*g 

where  E  is  given  by  the  eigenfunctions. 


3.3.43 

These  ex- 
the  fields 
.  It  fol- 

3.3.44 


IV .  Modifying  the  Expressions  to 
Account  for  Gain 


Chapter  IV  addresses  the  problem  of  generalizing  the 
previous  development  so  that  the  expressions  can  account 
for  the  presence  of  gain  in  the  resonator.  It  is  noted 
here  that  the  method  set  forth  here  is  not  the  only  way 
to  include  gain  in  mode  analysis  (Ref  8). 

In  this  thesis,  the  method  taken  to  include  gain  in 
the  preceding  development  will  vequire  two  changes  in  that 
development.  The  first  is  that  the  fundamental  cylindrical 
wave  is  no  longer  assumed  to  be  of  unit  amplitude.  In  the 
even  case,  to  which  consideration  will  be  limited,  g(x) 
is  then  assumed  to  be  of  the  form 

N 

g(x)  =  h  +  2  cnHn(x)  4.1.1 

n  =  l 

where  h  is  the  amplitude  of  the  basis  wave,  in  general  not 
equal  to  1,  which  will  be  determined  later.  The  second  mod¬ 
ification  made  is  to  include  a  gain  factor,  e^9(y)L  f  i n 
the  kernel  of  the  diffraction  integral.  The  integral  equa¬ 
tion  then  becomes 

1  —  - i t ( v-—) 2 

Xg(x)  =  /i  t/n  /*e2g(y^L  e  m  g  ( y )  dy 
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N 

X(h  +  X]  cnHn(  x)  ) 

n=  1 


/TE/7  />'*>  L  e'U(y^)2(h+£CrlHn(y))dy 
-1  n~l 


4.1.3 


The  best  value  of  h  can  be  found  from  a  relaxation  pro¬ 
cess  wherein  gains  are  assumed  to  be  equal  to  losses.  First 
however,  the  equations  require  a  round  value  as  a  starting 
point. 

To  determine  that  rough  value,  it  is  assumed  that 
there  is  a  uniform  intensity  across  the  laser  cavity,  in 
particular,  at  the  output  plane.  Thus,  the  gain  there  is 
affected  in  a  similarly  uniform  manner.  If  the  gain  med¬ 
ium  is  homogeneous,  then 


g(x)  = 


_2j> 


1  +  2 


nrr 

sat 


=  g 


4.1.4 


where  g0  is  the  small  signal  gain,  and  I  t  is  the  sat¬ 
uration  intensity,  both  determinable  from  actual  laser 
parameters.  I(x)  is  multiplied  by  2  since  there  are  in¬ 
tensity  contributions  from  two  waves,  one  propagating  in 
each  direction. 

If  the  uniform  intensity  across  the  feedback  mirror  is 
If  ,  then  the  feedback  power  is  given  by 


Pf  =  I f • 2ad 


4.1.5 
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4.1.13 


I(x) 
*  sat 


( 


ia  A  . 

1  n  m 


■) 


1. 

2 


s  q°i  .  i 

In  m  2 


4.1.14 


This  ratio  is  the  ratio  of  the  intensity  on  the  mirror  to 
the  saturation  intensity,  and  it  will  be  considered  as  the 
relative  intensity  of  the  fundamental  cylindrical  wave. 
Therefore , 


in  a  first  approximation.  This  will  give  a  rough  starting 
point  for  h  from  which. the  equations  can  begin. 

Cons ide ring  the  new  integral  equation,  4.1.3,  in 
light  of  the  first  stationary  phase  approximation,  it  is 
seen  that  a  term  has  been  added  to  the  various  q(y)'s 
From  the  approximation  it  is  then  concluded  from  the  re¬ 
sults  of  Appendix  A,  that 

N 

A(h+cnHn(x))  -  he2g^y°  ^L+he2g^  1^LH1  (x)  +  ^  e2g^y°  ^anFn+1(x) 

n  =  1 

N  _  b  N  _  N 

+  2e2g{yo)bnGn  +  i(x)+F1(x)2cnHn(l)e2g(1)L+G1(x)X;cnHn(-l) 
n=l  n=l  n=l 

e2g(-l)L  4.1.16 


44 


Since  the  intensity  profile  is  even,  in  this  case  it  is  as¬ 
sumed  that  the  gain  function  "g  is  even  too,  and  the  ap¬ 
proximation  then  simplifies  to 


A(h+EcnHn(x) ) 


he29/(y.)L+t,e2g(mHl(x) 


n=l  n=l 


+H1(x)e2g(1)L£cnHn(l) 
n  =  l 


4.1.17 


Equating  coefficients  of  Hn(x)  »  n  f  1  shows  that 


Acn+1  =  ane 


Zg(y°)LFn(x)+bneZg(y“)LGn(x) 


4.1.18 


However,  yf  and  yo  themselves  are  now  functions  of  x 
and  therefore,  the  sequential  arguments  leading  up  to  an 
eigenvalue  polynomial  can  no  longer  be  made. 

In  order  to  build  that  polynomial,  one  more  simplify¬ 
ing  assumption  is  made,  that  being  whatever  intensity  fluc¬ 
tuations  present  across  the  output  plane  exist,  their  ef¬ 
fect  on  the  gain  is  negligible.  The  gain  factor  is  then 
assumed  to  be  a  constant,  for  all  points  across  the  reson¬ 
ator. 

Defining  the  gain  factor 


where  g  is  given  by  4.2.4,  then,  the  equation  becomes 

N  N 

X(h+£cnHn(x))  =  h?  +  h5H1(x)  +  2^cnHn+l(x) 
n=l  n=l 

N 


+  H1(x)^CcnHn(l)  4.1.20 

n=  1 

Equating  coefficients  of  Hn(x)  now  gives 

cn+ 1  =  \  cn  4.1.21 

c»  >  f  c,  4.1.22 

cn+l  ■  "•1-22 


and  in  turn  it  is  seen  that 


c  -  r  — 

Cn'c 


c  (^)N 
CN  V 


4.1.23 


Equating  coefficients  of  H  (x)  yields  as  before 

N 

AC!  =  h?+£?cnHn(1)  4.1.24 

n  =  l 


Equating  constant  terms  shows  that 

Ah  =  h£+£cNHN+1  4.1.25 
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4.1.26 


AH-hg  = 
CHn+i  n 


M*rO. 

CHN+1 


c 


n 


Therefore,  from  4.1.27  and  4.1.23, 


U-Oh 

HN+1 


Substituting  this  into  4.1.24  yields 


X(i)N-'  H  ( 1 ) 

t.  e^N+i  n 


N+l 


n  =  l 


which  simplifies  to 


,  N  JL  ,  N-n 

l)(K)-iH  Hn(1) 


n  =  1 


and  the  polynomial  is  then  given  by 

AN(  x-c  )  =  ?N+1HN+1  +  (A-S)X1  ^N"n?N5n"NHn(l) 

n  =  l 


CN+1Hn+1+(X-5) 


„d) 


n  =  1 


4.1.27 


4.1.28 


4.1.29 


4.1.30 


4.1.31 


4.1.32 


The  roots  of  this  polynomial  can  be  found  using  the 
same  method  as  before,  which  will  then  be  rough  approximations 
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to  mode  eigenvalues.  If  the  model  is  to  be  correct  for  a 
steady  state  resonator,  gain  should  just  balance  loss  in 
that  resonator,  implying  that 

u'(x)  =  u(x)  4.1.33 


and  therefore 


Y  =  1  4.1.34 

If  y=l  ,  then  from  2.1.32 


X  =  /m  4.1.35 

Using  this  condition,  4.1.35,  then  h  can  be  modified 
until  the  lowest  loss  mode  has  an  eigenvalue  equal  to  /in  . 
The  value  of  h  that  allows  this  should  then  be  the  most 
reasonable  value  of  h 

When  h  is  found,  the  proper  gain  factor  is  in  turn 
found  by  substituting  h2  for  the  intensity  ratio  in  4.1.4 

To  extend  these  solutions  beyond  the  shadow  boundaries 
one  merely  has  to  multiply  the  constants  an  and  bn  in 
the  expressions  derived  in  the  last  section  of  the  previous 
chapter  by  £  ,  and  change  the  constant  factor  to  h 


8 


V.  Results 


Implementation  of  Code  and  Result  Check 

The  expressions  developed  in  chapters  three  and  four 
were  incorporated  into  a  CDC  Fortran  IV  program,  BARC,  which 
was  organized  into  two  basic  sections.  The  first  section 
included  development  of  the  coefficients  of  the  eigenvalue 
polynomials  3.2.15  and  3.2.26,  the  computation  of  the  roots 
through  IMSL  routine  ZCPOLY,  the  computation  of  the  weight¬ 
ing  constants  cn  according  to  3.2.16  and  3.1.23,  and  the 
preliminary  eigenvalue  expressions  based  on  3.1.2.  The 
second  part,  in  a  separate  subroutine,  implemented  the  ex¬ 
pressions  developed  in  the  third  section  of  chapter  three: 
the  eigenvalue  expressions  valid  for  all  x  .  The  program 
was  then  run  for  various  cavity  parameters  and  the  results 
were  compared  with  the  results  of  other  programs  (Ref  6:1536; 
Ref  8:239). 

Table  1  represents  a  comparison  of  eigenvalue  moduli 
resulting  from  the  program  developed  in  this  work,  and 
those  from  the  Moore  and  McCarthy  program. 

These  results  are  for  a  cavity  with  magnification  of 
2.9  and  an  effective  Fresnel  number  of  16.4  .  The  solu¬ 

tion  compared  is  that  of  the  even  parity  case. 

It  is  seen  that  the  two  codes  predict  modes  with  very 


TABLE  1 

Mode 

Mod  BARC 

Mod  m&mc 

1 

1.040105 

1.040105 

2 

.  .625501 

.625501 

3 

.606668 

.606668 

4 

.496561 

.496561 

5 

.467285 

.467285 

6 

.  157664 

.  139999 

s imi 1  a  r  1 osses  ,  s i nee 

loss 


5.1.1 


and  the  A*A  values  are  all  very  close. 

Figures  3  through  8  are  included  to  show  results  of 
eigenfunction  intensity  plots  over  similar  ranges  for  the 
Moore  and  McCarthy  program  and  program  BARC.  Figures  10 
through  15  show  comparison  between  BARC's  results  and  those 
published  in  reference  6  (Ref  6:1536-1539).  In  both  cases, 
through  visual  comparison,  program  BARC  produces  results 
that  are  very  similar  to  results  from  previous  methods. 

This  indicates  that  BARC  produces  valid  re suits,  at  least  to 
the  extent  that  the  previous  methods  are  valid. 
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Bare  Cavity  and  Gain  Results 


After  the  validity  of  the  code  was  ascertained,  the 
code  was  modified  according  to  the  expressions  generated 
in  chapter  four.  Figures  14  through  21  illustrate  results 
obtained  for  a  bare  resonator  of  magnification  2.9  and 
equivalent  fresnel  numbers  of  15.863  and  16.4  .  These 
parameters  are  chosen  to  facilitate  mode  separation  com¬ 
parisons  later  in  this  section.  Figures  22  through  29  il¬ 
lustrate  results  of  a  loaded  cavity  of  the  same  configura¬ 
tional  parameters  but  containing  a  gain  medium  of  small 
signal  gain  5%cm_1  and  cavity  length  of  200cm.  This 
group  of  plots  allows  comparison  between  bare  and  loaded 
cavity  cases.  It  is  seen  that  this  particular  resonator 
model  predicts  that  loaded  cavity  modes  have  nearly  the  same 
intensity  profiles  as  bare  cavity  modes,  differing  only  by 
a  scale  factor. 

At  first  glance  this  seems  reasonable,  since  in  the 
bare  cavity  case,  the  whole  eigenfunction  was  based  on  a 
wave  of  unit  relative  amplitude,  and  slight  modifications 
on  that  wave  by  diffraction  supplied  by  the  oscillatory 
functions  Hn(x)  •  In  the  loaded  cavity,  the  eigenfunction 
is  also  based  on  a  wave  modified  by  the  same  functions,  only 
the  relative  amplitude  of  that  wave  is  no  longer  unity. 

Thus  it  seems  likely  that  the  profile  would  look  moderately 
similar  in  both  cases. 
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CONSTANT 


In  closer  analysis,  the  mode  for  the  bare  resonator 


has  eigenfunctions  given  by  this  expression 


9b(x)  = 


5.2.1 


while  the  loaded  case  has  eigenfunctions  given  by 


gL(x)  =  h  +  XXLHn(x) 


5.2.2 


Similarly,  the  expressions  for  the  weighting  constants  are 


^  b  { X- 1 )  ,N-n 


5.2.3 


C  L  =  <X-g)  h(^)N’n  i 


5.2.4 


the  eigenvalue  polynomial  in  the  loaded  case  is  (4.1.32) 


XN(X-C)  =  gN  +  1HN+1  +  (X-?)ExN-nCnHn(l)  5.2.5 


Dividing  through  by  £  yields 
N  N 

(r)  (r-D  =  Hn+1  +  (|-1)£  XN_nHn(l)  5.2.6 


which  becomes  identical  to  the  bare  cavity  polynomial  (3.2.15) 
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as  £-*-1  .  This  condition  will  be  fulfilled  when  h  ,  the 

intensity  ratio,  becomes  very  large  as  seen  from  4.1.19  and 
4.1.4  .  In  turn  H  is  then  seen  that,  as  £->-1 

CL  -  C  b  5.2.7 

n  n 

and 

gL(x)  -*  hgb ( x )  5.2.8 

From  this  it  is  concluded  that  in  the  well  saturated  case, 

or  when  the  ratio  of  the  actual  intensity  to  the  saturation 

intensity  is  much  more  than  one,  the  field  distributions 

and  hence  the  intensity  profiles  will  equal  those  of  the 

2 

bare  cavity  case  multiplied  by  h  and  h  respectively. 

Tables  2  and  3  are  presented  to  illustrate  and  compare 
mode  separation  properties  of  a  loaded  and  a  bare  cavity 
for  three  different  equivalent  Fresnel  numbers.  The  para¬ 
meters  chosen  were  a  a  magnification  of  2.9  and  N^'s  of 
16.874,  16.4  and  15.863.  These  were  shown  in  reference  6 
(Ref  6:1534)  to  be  points  of  least,  greatest  and  then  least 
loss  and  next  to  lowest  loss  eigenvalue  moduli.  It  was 
thought  that  since  the  lowest  loss  mode  eigenvalue  was  forc¬ 
ed  to  the  same  constant  value  at  each  Fresnel  number,  negat¬ 
ing  any  quas i peri odi ci ty ,  the  higher  loss  modes  might  also 
lose  quasi  periodicity.  The  numbers  presented  show  that  the 
higher  loss  modes  do  maintain  their  quasi  periodicity. 


TABLE  2 


BARE  RESONATOR 

Mod 

(X) 

Mode 

Nf=15.863 

Nf =16 . 400 

Nf=16.874 

1 

0.8543652 

1.040102 

0.8922496 

2 

0.8508141 

0.6255715 

0.7785354 

3 

0.5385818 

0.6067205 

0.5400256 

4 

0.5049350 

0.4966156 

0.5290538 

5 

0.4737932 

0.4673182 

0.4752758 

6 

0.1718837 

0.1646309 

0.1593562 

TABLE  3 

LOADED  RESONATOR 

Mod 

(X) 

Mode 

Nf= 15 . 863 

Nf=16.4 

Nf=16.874 

1. 

1.702922 

1.702965 

1.70294 

2 

1.695844 

1.024215 

1.485906 

3 

1.073502 

0.9933517 

1.030688 

4 

1.006437 

0.8130828 

1.009747 

5 

0.8718129 

0.7651157 

0.9071073 

6 

0.3425989 

0.2695415 

0.3041459 

The 

h's  required 

to  adjust  X  to 

^/m  were  2.6899 

3.1051,  and  2.5979  for  N^=15.863,  16.4,  a  rd  16.874  respectively 
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V I .  Concl us  ion  a nd  Recommendations 


Concl us i on 

The  primary  conclusion  of  this  thesis  is  that  pro¬ 
gram  BARC ,  written  according  to  expressions  developed 
along  Horwit'z  analysis,  produces  valid  results.  The  pro¬ 
gram  allows  analysis  of  even  and  odd  parity  mode  solutions, 
more  general  than  Moore  and  McCarthy's  program,  and  also 
allows  field  calculation  beyond  the  shadow  boundary. 

Incorporation  of  gain  considerations  into  the  program 
to  allow  analysis  of  a  loaded  strip  resonator  has  been  done. 
After  modification  the  program  produces  results  from  which 
a  second  conclusion  can  be  drawn,  that  being,  for  this  par¬ 
ticular  model,  mode  intensity  profiles  in  a  loaded  strip 
resonator  are  essentially  the  same  as  those  predicted  for  a 
bare  strip  resonator.  It  is  also  concluded  that  mode  losses 
as  function  of  equivalent  fresnel  number  continue  to  exhibit 
quasi  periodicity  in  the  loaded  case. 

Recommendati ons 

The  computer  program,  as  it  stands,  predicts  some  very 
basic  results  about  modes  in  an  unstable  resonator.  There 
is  no  doubt  that  the  scope  of  the  program  and  the  model  upon 
which  it  is  based  can  be  broadened  considerably.  As  it 
stands,  it  could  be  used  to  examine  a  full,  or  more  complete 
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range  of  resonator  parameters  either  loaded  or  bare. 

The  program  should  be  used  to  explore  mode  separation 
in  the  loaded  cavity  case.  Mode  separation  could  be  ex¬ 
amined  for  a  range  of  Fresnel  numbers  as  has  been  done  for 
bare  cavities.  (Ref  6) 

The  model  given  here  could  be  modified  to  account  for 
a  non  uniform  gain  function.  To  do  this  a  new  series  of 
Hp's  might  be  developed  through  asymptotic  analysis  of  the 
gain-modified  kernel.  Another  method  might  be  the  use  of 
matrix  methods  to  solve  the  eigenvalue  equation. 

This  existing  method  could  be  applied  to  resonators 
with  circular  mirrors  per  Ref  2,  and  gain  then  included 
in  that  case. 
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Appendix  A 


This  appendix  will  employ  the  stationary  phase  approx¬ 
imation  to  simplify  the  resonator  integral  equation  into  a 
workable  expression.  The  derivation  starts  with  the  func¬ 
tions. 


and 


F(x,t) 


-1  e-it(l-x)2 
2/TirT  1  -  x 


A1 


G(x,t) 


-1  e"1 1+x ) 2 

2/TiT  1+x 


A2 


These  are  modified  by  letting 


and 


£ 

where 


F„(x) 


=  f(£tt  » 


n- 1 


Gn(X) 


n 


K  =  0 


A3 


A4 


A5 
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and 


m  =  magnification 
i  =  /^T 

Thus  it  is  seen  that 


F, (x)  = 


-1  e 


-  i  t  ( 1  -  x  /  m ) : 


2/vttT 


1-x/m 


A6 


F,(x)  = 


./TTTTiF  e-it(1'x/m2)  /1+1/m 


F_  (  x  )  = 


and  similarly 


Gn(x) 


2  i  t 

1-x/m2 

that 

- /m  , 
■V  n-  1 

-it(l-x/mn)2/mn-1 

e 

2/TTrT 

l-x/mn 

-VvT 

- i  t(  l  +  x/ml1)2/mn_1 
e 

2/lTTt 

l  +  x/m11 

A7 


A8 


A9 


Now,  it  is  recalled  that  the  working  form  of  the  inte¬ 
gral  equation  is 

1 


U(x> 

-1 


dy 


A10 
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In  the  even  part 


N 

'  f(x)  =  1+  £  <anFn(x)+bnGo{x))  A11 

n=l 

ATI  is  substituted  into  the  integral  equation,  which  now 
becomes 

X(l+E{anFn(x)+bnGn(x)})  = 

it  1 

—  j  elt(y-x/m)2(l+I{anFn(y)+bnGn(y)})  dy  A12 

-1 

When  expanded  once,  the  right  side  becomes 

=  ^pE  j  e-it(y-x/m)Zd y+7?  2  /  e-it(y"x/m)2 

-1  n=l  -1 


The  first  term  is  called  I0  .  Then 

I,  *  VF  /  e-U^-x/m)!dy 

-1 


A13 


A14 


This  is  now  considered  in  light  of  the  first  order  approxi¬ 
mation  to  the  method  of  stationary  phase  which  states  that  if 


A 1 5 


I  =  /  q(y )e 

a 


■1  tp(y) 


dy 


then 


I  a  e_i7r/4a(vn)eitp(yo)  f  -2f 

q[y°,e  "l/tp  (y« ) 


+  r 


Sib)  e-itp{b)  _  g(a)  -itp(a) 
P  (b)  p'(a ) 


Where  y0  is  such  that 


A16 


P  ^  ( y  o )  =  o 


A17 


It  is  seen  that 


q(y)  =  1 


p(y)  =  (y-x/m)2 


p'(y)  =  2 (y-x/m) 


p~(y)  -  2 


A17.1 

A18 

A19 

A20 


And  that 


yo  =  x/m 


A21 
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] 


The  first  term  in  the  sum  of  integrals  is  now  called  I 
1 

Ij  =  STTp*  J  e-lt(y-x/m)2  (ajFjty)  +  bjGjfy))  dy  A30 

-1 


Explicitly,  this  becomes 


1 

/ 

-1 


r  |  -a,  -i t( 1-y/m) 2 

Ij  =  /TT77  J 


2/TiTt  1-y/m  2/TriT 


-i t( 1+y/m) 


1+y/m 


•  i t (y-x/m) : 


A31 


Upon  separation  the  result  is 

1 


I  -  /TT7?  f  — i-  e 
J  2/TTt 


-a,  -it| (l-y/m)2+(y-x/m)2 


-1 


1-y/m 


dy 


+  /i  t/TT  J 


"bl  e-it| ( l+y/m)2+(y-x/m)2 


-1 


2/Tit 


1+y/m 


A32 


Upon  consideration  of  the  first  part  of  this,  it  is 
seen  that 

"  T^7i  ft33 


p ( y )  =  (l-y/m)2+(y-x/m): 


A34 


p'(y)  = 


-  (I-1)  +  2 ( y -— ) 
m  m'  m' 


A35 


p"(y) 


+  2 
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A36 
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P(y>  -  (^‘♦(y-i)* 


A42 


P'(y)  *  t(l+^)+2(y-i) 


m' 


A43 


p"(y)  -  *  2 


A44 


and  solving  for  y0  ,  the  result  is 


0=1+  +  2y  o  -  — 

m  mz  J  0  m 


2x 


A45 


y  o  ( — *r+ 1 )  -  -  ~r 

J  m  m  m 


A46 


y.  ■  <=  -  =>■  1 


m  m'  ,.T 


A47 


Substitution  of  these  expressions  into  the  second  part  of 


A32  results  in 
r 


-b 


1 


-iir  -i  t 
4 


2/ifft 


/i~t/ff<e 


tt  e 


"V^  m  V  (l4r  ”  m  ”/ 


tl+^2“ 
m* 


i  m  i.l  mm 

1  + — y 

nr 


( 1+-)  *  {— ( 1+— }  +  2  ( 1  -— ) )  (1-— )*(— (1-— )+2(-l~— )) 
v  m'  vmv  nr  '  m  v  nr  'nr  nr  v  nr' 


>  A48 
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Considering  the  denominator  in  (A-40)'s  stationary  phase 
point  contribution,  it  is  seen  that 


1- 


l 

m 


4r(? 


1+rrr 

m 


+ 


1- 


x+1 

(m+^-) 


A49 


And  similarly  in  (A-48)'s  stat  phase  point  cont,  it  is  seen 
that 


1+ 


l 

m 


1 


(- 

m 


1- 


x-1 

(m+i)  »m 


A50 


Also,  it  is  seen  that  the  argument  of  the  exponent  in  A-40's 
stationary  phase  contribution  can  be  simplified  as  follows 


A51 


A52 


A53 


A54 


A55 
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2 


(l+£r> 


m  +  m  '  m 


<1+5F> 


(^r+l) 

'  m  *■ 

7~ 


(l+^O 

m 


A65 


(jjT+l)(l+jr) 

<1+St) 


I1-  p-> 

1+-L- 

m 


A66 


Substituting  these  simplified  expressions  into  A-32, 
find  that 


we 


I ,  =  -/ft/ rr 


1 


2/TFt 


-i-rr/4 


IT 


-it(l-£r)2/l+^ 
m  m 


1  - 


x+1 


"("V 


ble 


-  i  t  ( l+£) 2  /  1+^-2-v 


x  - 1 


/  i  t  /  n 


1  i 


m(m+rn) 


2/TfTt  t 


-it((l-I)%(l-l)2)  -it((l+I)2+(-l-£)2) 

aje  a,e 


t 


( ^  >  <  1+i>  (ir(  1+i'+2<  -1-! 

U(l+I)!  +  (l-i)2  _  -U(l-i)2+(-l-i)2  ] 


)) 


ble 


ble 


^  A67 


Now  the  first  part  of  this  expression  can  be  simplified  as 

fol 1 ows : 


1- 


1 


x+1 

m(m+4) 

m 


1- 


1 

x+1 

ITT  +  1 


A68 
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1 

nr  + 1  -  x  -  1 
m 1  + 1 


A69 


m2  + 1 
m2  -  x 


A70 


A71 


And  the  second  part  can  also  be  simplified 


1  =  1 

i  ■  x- 1  m z  + 1  +  x -T 

nF+T  m'2  +  l 


A72 


rrF  +  1 
in  2  +  x 


1+ 


1+ 


A73 


A74 


If  these  simplified  expressions  are  substituted  into  A67,  the 
result  is 


I ,  a  -  /  i  t  /  tt  — - —  e  " 77 !  ^ 
1  2/i  TTt 


ale 


•  i  t(  l-^-)2/l+ir 
nr '  m 


1-^ 


m 


ble 


n+Li 


1+^ 


1-  /i  t/ tt 


I  1 


2/TiTt 


2 1 
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ale 


ala 


-it((lU)2  +  (-l-|;)!) 


(l.Ijf.I+lr+l.l) 

'  m' x  m  m  m 

-u((i4)zmi-S)!) 


bie 


m 


ble 


(1+I)(I+1T.+  i_2L) 

'  nr  'm  m  nr 


(1_I)  (I-i,--  l-x) 

'  nr  vm  m  m' 


A75 


After  a  few  sign  manipulations  and  cancellations,  it  is 
seen  that 


i  _ / a,e 

I,  =  — —  /T+T7'm7'  ( 


1  2/rfrt 


1  - 


-i  t(  l+^z-)2/l  + 


ble 


m 


1  + 


1 


4  1  7T  t 


a  e’it(1"m)  -it(l-£)  'it(1+m)2  -it(-l- 

1  e  !L_  +  !il  e 


( 1+i)  ( 1+1+-4+-) 
v  nr  v  m  nr  m' 


bxe 


n<ui>!  -itd-i) 
m 


(1+I)(1+I+lT.ii) 

nr  m  nr  m 


„  -U(-l-i)2 

b, e  '  m' 

1 _ e  _ _ 

(l-Ijd-ifly+i) 
v  nr s  m  nr  nr 


r  A76 


The  stationary  phase  point  contributions  are  seen 
immediately  to  be  equal  to 


ajF^xJ+bjG^x)  A77 

In  the  end  point  contributions,  the  denominators  must 
be  approximated  and  terms  of  ^  or  higher  order  be  neglect¬ 
ed.  When  done,  the  end  point  contributions  appear  as 


9 


3  lx  3  I*— 


Il  -  a1F2(x)+b1G2(x) 

+a1(Fl(x)F1(l)+G1(x)F1(-l)) 

+b1(F1(x)G1(l)+G1(x)G1(-l))  A81 

The  second  term  in  the  sum  of  integrals  is  now  con¬ 
sidered  as  I2  : 


98 


/i  t /  tt  f  e 

-1 


-it(y-£) 

m  (a2F2(y)+b2G2(y))  dy 


1  X  2 


=  /T t/n  J  e 


-1 


-a,  li  -ItU-Ja-l'/Hi, 
£.  m  e 


2/i-rrt 


*  -i*- 


h  1  +  1  -it(l+*4)2/l+m- 

b2  1 V  e  m  m 


2/TirE 


1  1 


1 


dy 


A82 


Considering  the  first  term  of  this  exp>ession,  it  is  seen 
that 


q(y) 


1 


1-Jr 

m 


A83 


p(y)  =  (y-— )  2  + — r — ( 

w  m'  j  +  l  m 

m7” 


A84 


p'«»>  ■  2(^S)-37A-T!1-^) 


“  "2<1+iF> 


A85 


p"'(y)  =  2+ 


mM1V) 


A86 


And  solving  for  y0  it  is  seen  that 


0  =  y  o 


(m2(i+^r)) 


A87 


99 


Considering  the  second  part,  it  is  seen  that 


q(y)  = 


A92 


P(y>  -  <y-^)°  *  -V  d*t) 


l^~T 
m 


A93 


p-(y)  =  2(y-i)  +  2  1  (1+L, 


1+t^T 

m 


A94 


P~(y)  =  2+ 


^(1^) 


A95 


Solving  for  y 


0  J 


0  =  y  -i+ij.  -L-d+ii) 

"  m  ir  ,,1  v  nr' 


1+- 


nr 


A96 


0  =  y0--+-T-  — —  +  - ^~r- 

J  m  nr  uL  -4/  ,.l 


l+£r  mMl+V) 
nr  '  m  ' 


A97 


y0  fl  + 


(  mMl+jU*)) 


m 


1 


m2  ( 1+m7’) 


A98 


and  finally 


=  fu  - 


m2{1+ih‘)A1  +  — 7~T 


m'*(1Vr)' 


A99 
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-it((y,-i)Wl+£f)Vl4r) 


b2e  e 


1  1  a2e 


‘♦J*  J 

m7" 


2/TT?  L  (l-±r)(2(l-±) - S— d-^-) 

m  m  m2(l+ir) 

m 

1+ — n" 

m  *■ 


4r)  ( 2  ( -  l-£) - - — ? ( 1+K) ) 

m 

-i  t(  ( 1-— )  2+ — r-d+^r)) 
v '  nr  1  +  1  m 

rri7" 

lr<a‘1-s)+ir-  7Ti-<1+ir>> 


1+iF 

•> 

1  m* 

-)(2(-l-%)  +  ?-  (1--U-) ) 

m  rn2(l+^r)  m 


rrnr  /i  t / tt 


A10 1 


The  denominators  of  the  stationary  phase  point  contribution 
terms  must  be  simplified: 

In  the  first  term, 


A102 
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And  in  the  second  term, 


Now,  considering  the  argument  for  the  exponent  in  the  first 
part,  it  is  seen  that 


=  1  (l  .  X  \/  1+1/m2  y 

1+p-  \  ^/V  +  l/m^l/mV 


the  total  argument  is 


1 


+  ( 1  + 


A123 


A124 


A 125 


A126 


A127 


second  half. 


A128 
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109 


Thus  the  stationary  phase  point  contribution  simplifies  to 


'  ^  -It(d-Sr) 


- /T  /l+l/m^+l/mL 
2/i  it  t 


i 


a^e  e 


1  - 


A 145 


- 1  it 


And  this  is  seen  to  be,  since  e  /7  =  , 


a2F3(x)+b2G3(x) 


A146 


Now,  considering  the  end  point  contributions,  the  denomin¬ 
ators  are  expanded  and  the  2's  are  factored  out  to  get 


X\  2 


-/1+1/m 

4/TFt 


U/  V\ 


a2e 


it/-,  1  i2 


(l-  -4)d-  -  -  — — (i-  ^)) 

^  m  m2+l  m 


•itd+i) 


U  (1+^-)2 


a2e 


-ef<1+s->2 

b2e  e  ^ 


b2e 


-  d“(1-  ^)2  ^ 

in7" 
e 


+  iFTT(1'SPr,) 


A147 


Once  more  in  the  denominators,  terms  of  1/m  or  higher  are 

111 


neglected,  and  it  is  seen  that  the  end  point  contributions 


- /TtTrr  7 

4/Ttt  t 


-it(i-j)2  -  -4-{i4)j 

m  t  ,  1  m 


i  J  a2e 


*  l  c •  =0  e-o 


1  m7"  .  in7" 

a9e  e  b?e  e 

+  — -  +  — - - - - - : - 


M)  W 


0*0  (•♦* 


a2F1(x)F2(l)+a2G1(x)F2(-l) 

-b2F1(x)G2(l)+b2G1(x)G2(-l)  A150 

Thus  it  is  seen  that,  after  adding, 

1 1  +  1 2  =  a1F2(x)+b1G2(x)+a2F3(x)+b2G3(x) 

+  F1(x)(a1F1(l)+b1G1(l)  +  a2F2( l)+b2G2( 1) ) 

+  G1(x)(a1F1(-l)+b1G1(-l)  +  a2F2(-l)+b2G2( -1) ) 
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A151 


From  this  it  is  concluded  that 


N  X.  -it(y-i)2 


f\  t/ir 


e/ 

n=l  -1 


anFn(J')+bnGn(J')) 


N 

> 

n  =  l 


■  EtanWX>  + 


tFl<x>  E  <anFn(1)+,‘»Gn<,» 

n  =  1 
N 

tGl<x>  E  UnFn<-1)+bnGn<-1>> 
n  =  l 


And  in  turn,  adding  the  I0  term,  it  is  had  that 

X) 

(l+E<anFn(y)+bn(FnCy)J) dy 


1_ 


/i  t/  it  /  e 
-1 


=  1  +  FjU)  +  Gj(x)  +  2  {anFn  +  l(x)  +  bn( 

n  =  l 


n  +  1 


N 

> 

n  =  1 


F  GjU)  ^(anFn(-l)  +  bnGn(-l)} 


A152 


(x)  > 


A153 
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Appendix  B 

This  appendix  simplifies  the  definite  integral  pro¬ 
duced  by  assuming  in  2.1.16  that  a£  is  effectively  in¬ 
finite. 


oo 


~2L [gixi2  +92x22  -2xlx2]  "2T  [9lxlZ  +  g2x2'2xlx2] 


ik 


dx: 


,  /•  •k[9lxl2  +  92x22-2xlx29lxl2  +  92x2-2xlx2] 

=  ir  J  e  dx; 


B 1 


CO 

-ir/- 


'2L  [291xi2"2(x2  +  x2^xl  +  92^x22  +  xp] 


dx: 


B2 


"2L  92 ( x22  +XI ^  r  " 1 L [9 1X 1 2 ” ^ X2  +  X2 ^ x l] 

xre  J* 


dx: 


B3 


The  square  is  completed  in  the  exponent  by  adding  and 
subtracting  b2/4a  ,  or 


ik. 


XL 


/' 


“  2L  9  2  (X2Z+X2  ^  °r  '  L  9lxl2"^x2  +  x2 


(x«+Xp)  "] 


ik  ( x2  +  x2 ) 


L  4g . 


dx  J 


B4 
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115 


then 


=  Yj-  ( 1+i )  A/2  /L/-g  jK  e 


B 14 

B15 

B 16 


which  is  the  final  kernel. 


% 
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IP 


Appendix  C 


This  appendix  simplifies  the  integral  equation  in 
2.1.29  into  the  final  form. 


1  -  i  tt  F  [  g  (x2+y2)-2xy] 

yu(x)  =  /  /TTu(y)e 

-1 


Cl 


Let 


u  ( x )  =  g  { x )  e 


-k(m  -1)x 


C2 


since 


«-irF 


C3 


then 


yg(x)e 


=  /TT  /g(y)e 

-1 


1  -i7rF[g(x2+y2  )-2xy] 


However,  since 


/g+l  -  /g- 1 


dy 

C4 

C5 


-  (/qTT  +  /g^T)2 
g+l-g+1 


C6 


117 


2g  +  2  / (  q'+ 1 )  (  q -  1 ) 
2 


C7 


g  +  /g^T 


C8 


so 


m2  =  g2+2g/g2-l+g2-l 


C9 


m2  +  l  =  2g2  +  2g/g2 -  1 


CIO 


■  9!+g/g^T 


Cll 


If  both  sides  are  divided  by  m  ,  or  rather  one  by 
m  and  the  other  by  its  equivalent,  g+/g2 - 1  ,  the  re¬ 

sult  is 

m2  +  l  _ 

2'in  ”  9 


C 12 


this  can  then  be  substituted  into  the  integral  in  place  of  g: 


1 

yg(x)  =  /TV  I  g (y ) e 
-1 


■  i  7T  F  ■ 


m2  +  l  ,2j_ni2  +  1..2  -1..2  m2 


2m 


x2+^y2-2xy+^y2-i 


2m  J  2m 


dy  C 13 


i  _ilTF  fm-  +  my2-z*y\ 

/TT j  g (y ) e  dy 


C 14 


=  1 


f, 


I.  -luFfy-*)’ 


=  /TF  /  g(y)e 


dy 


C15 
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PROGRAM  BARC(  If! PUT, OUTPUT ,TAPK3= OUTPUT) 

REAL  NEQ ,  MAG ,  MS'JE!.’.'  (  5 1 )  ,  HSUPN!  5  1  )  ,  I  "TEN  L1  1000) 

COMPLEX  EYE  ,  COEF (51)  ,  Af! (51)  ,  LAND.;/.!  5  1) 

COMPLEX  CL(  5  1  )  ,  CONST!  51),  CDUM ,  A 1  ,  AN 2  ,  RTS  YE 
COMPLEX  FIELDXi  1000)  ,SIG,3!.'1  , ,  ROOT 
DIMENSION  LABEL (17)  ,  5TCF.EX ( 1  COC )  , sINDEX ( 51 )  ,  PL0C0N(51) 
DIMENSION  FUKVALC. -0,410)  ,  PLOF'J!.(  BO) 

DATA  LABEL/ 1 7( 1 CH  )/ 

C 

C 

C  THIS  PROGRAM  COMPUTES  RESONATOR  MODE  EIGENVALUES  AND 
C  SUBSEQUENT  EVALUATES  INTENSITY  VALUES  FOP.  POINTS 
C  ACROSS  THE  OUTPUT  PLANE  OF  A  STRIP  LASER  RESONATOR, 

C  THE  PROGRAM  DEALS  WITH  EITHER  A  BARE  OK  LOADED 
C  CAVITY,  USER'S  PREFERENCE . 

C  -OUTPUT  CONSISTS  OF  AN  EIGENVALUE  LIST,  WITH  PHASE 

C  AND  MAGNITUDE,  FIELD  VALUES  FOR  A  SELECTED  NODE 

C  EITHER  ON  OR  OFF  THE  MIRROR,  PLOTS  OF  FIELD  SERIES 

C  FUNCTIONS  OR  WEIGHTING  CONSTANTS,  AND  PLOTS  OF  INTENSITY 

C  ACROSS  THE  OUTPUT  PLANE  WITH  EITHER  LIMITED  OR  EXTENTED 
C  RANCE. 

C  COMPILED  CODE  NEEDED  AROUND  110000  OCTAL  TO  LOAD. 

C 

C  INPUT  QUANTITIES  ARE  AS  FOLLOWS: 

C 

C  MAG  =  CAVITY  MAGNIFICATION 

C  NEQ  =  EQUIVALENT  FRESNEL  NUMBER 

C  MTEST1  =  FIELD  SOLUTION  PARITY  DESIGNATOR 

C  NBIG  -  DESIRED  It  TERMS  IN  FIELD  SERIES 

C  CAVLEN  =  CAVITY  LENGTH  IN  LENGTH  UNITS  FOR  LOADED  CASE 

C  GNAWT  r  SMALL  SIGNAL  GAIN  IN  PER  LENGTH 

C  H  s  AVERAGE  CAVITY  INTENSITY,  OR  EIGENVALUE  FORCING 

C  PARAMETER 

C 

C 

C  TO  TERMINATE  PROGRAM,  INPUT  MAG  =  0  OR  LESS. 

C  NOTE:  EVMAG  DENOTES  EIGENVALUE , MAGNITUDE  AND  EVPH 
C  DENOTES  EIGENVALUE,  PHASE. 

C 

C  THIS  PROGRAM  ALSO  REQUIRES  IMSL  ROUTINE  ZCPOLY  AND  PLOT 
C  LIBRARY  CCPLOT56X.  FINAL  COPY,  20  OCT  1980.  J  E  ROWLEY 

C 

c 

c 

99«  F0RMATCG10.3) 

LABEL( 1 ) = 10H  HEO= 

LABEL! 3 )  =  1 0H  MAG  = 

777  WRITE! 8,999) 

999  FORMAT! 1111 , IX, » INPUT  MAG,  NEQ,  AND  PARITY:  »,/) 

READ  * , HAG , NEQ , MTEST1 
IF(MAG.  l.E.O. )  GO  TO  383 
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WRITEC  8,88)  MAG  ,  NEQ 
IFCKTEST 1 . EQ . 0)  GO  TO  8 
WRITEC  8 , 977 ) 

GO  TO  9 

8  WRITEC  8 , 976) 

9  CONTINUE 

977  FORMAT (IX, *  PARITY  IS  ODD.  *,/) 

976  FORMATC IX, *  PARITY  IS  EVEN.  *,/) 

LABEL! 5 )  =  1 0H  NOPE  E 

LABELC6 )= 10HIGEN VALUE : 

C 

C  MSUPN(I)  =  MAG“(I-1  ) 

C  MSUBN(I)=1+1/ MAG*  *2  +  ...  +  1/KAG“C2*I-2) 

C 

MSUBNC 1 )  =  1 . 0 
MSUPNC1)=1  .0 
DO  10  1=2,51 
MSU  PN ( I )  =  M AG  * MSUPN ( I - 1 ) 

MSUBNC I )=MSUBN( 1-1 ) +1 ,/MSUPHC I) **2 

10  CONTINUE 
C 

PI=2 . * ASINC 1.0) 

EYE=CMPLXCO. ,  1  .) 

RTEYE=CMPLX( 1 . , 1 . )/SQRT(2.) 

DUM1  =  2* PI* MAG*  * 7  /  (  ”AG“2  - 1  . ) 

RNBIGr ALOGC  250*  NEC ) / ALOG ( MAG ) 

IF( RNBIG. LE . 50 . )  GO  TO  15 
WRITEC  6 , 998) 

GO  TO  777 

15  WRITEC 8, 996) RNBIG 

996  FORMATC IX, ‘CALCULATED  NBIC  =  * , G1 4 . 7 , ‘INPUT  INTEGER  CHOICE:* 
READ  * , N3IG 
WRIT E(8, 979) N BIG 
WRITEC  8 , 975 ) 

975  FORMATC IX, ‘TYPE  1  FOR  GAIN  CONSIDERATION  :  »,/) 

READ  * , IGAINQ 
WRITEC 8, 979 )IGAINQ 
IFC IGAINQ. NE. 1 )  GO  TO  5 
WRITEC  8 ,974 ) 

97*r  F0RHATC1X, ‘INPUT  LENGTH  AND  S-S-GAIN  IN  COMMON  UNITS  :  »,/) 
READ  * , CAVLEN , GNAWT 
WRITEC  8,97 1)CAVLEN,GNAV.'T 

C 

C  DIVIDE  INPUT  INTENSITY  GAIN  BY  TWO  TO  MAKE  IT  THE  FIELD  GAIN, 

C  WHICH  IS  WHAT  THIS  PROGRAM  ACTUALLY  REQUIRES 

C 

c 

GNAWTrGNAWT/ 2 . 

971  FORMATC IX, ‘INPUT  VALUES  ARE  :  *,2G1«.7,/) 

H = SQ  R  T  ( G  N  A  C T  *  C  A  V  L  E  N  /  A  L  0  G  C 1 1 A  G )  - .  5 ) 

WRITE(8,973)H 
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973  FORMAT ( IX, *H =*,014. 7, “INPUT  MODIFIED  VALUE  OR  0  TO  CONT  :  «,/) 
READ  * , HV  AL 
WRITE! 8,97?) HVAL 

972  FORMAT! IX, ‘INPUT  VALUE  IS  :  «,G14.7,/) 

IF! HVAL .  ME . 0 . )  Hi HVAL 

GAMMA=EXP(2*CAVLEN*(GNAWT/(  1  .+2*H**2))) 

GO  TO  6 

5  H=1.  $  GAMMAil. 

6  CONTINUE 
WRITE!8,993) 

993  FORMAT! IX, ‘INPUT  ZERO  TO  LIST  EIGENVALUES  :«,/) 

READ  * , LTEST 
WRITE! 8 ,979 ) LTEST 
LABEL! 1 3 )  =  1 0H  EVEN 
NDEG=NBIG+1 

IF! MTEST1 . EQ . 0 )  GO  TO  1  6 
LABEL! 1 3)  =  1 0H  ODD 

NDEG=NBIG 
16  T=DUM1 *  NEQ 

LABEL! 1 4) = 10H  PARITY 
C 

C  COMPUTE  COEFFICIENTS  OF  THE  POLYNOMIAL 
C  P(Z)=COEF(1 )*Z**NDEG  +  COEF !2) *Z** ! NDEG-1 )  +  ...  + 

C  COEF! NDEG) *Z  +  COEF! NDEG+1 ) 

C 

COEF! 1)=CMPLX!1.,0.) 

NCOEF=NDEG+1 
DO  25  1=1, NDEG 

AN1=RTEYE*2*SQRT(PI*T/MS'JBN(I) ) 

AN2r-T*EYE/MSU3.V!I) 

AN3=1 .-1 ./MSUPN! 1+1 ) 

AN4  =  1 ,+1 ,/MSUPN! 1+1 ) 

AN!I)  =  !CEXP(AN2*A.';3*52)/AN3  +  CEXP!  AN2*AN4«*2)/AN4)/AN1 
IF! MTEST1 . EQ. 0 )  GO  TO  25 

AN! I)  = ! CEXP!  AN2*  A.‘,'3 *  *2 ) /  AN3-CEXP!  AN2* AN4 **2 )/AN4  )/ AN  1 

25  CONTINUE 
IF!MTEST1.EQ.1)  GO  TO  27 
COEF! 2) =( AN! 1 ) -1 . ) ‘GAMMA 
DO  26  1=3 , NCOEF 

26  coef!  i )  =ga;:;:a*»  ( i-  i )  *!  an!  i-  i  )  -ahi  1-2 ) ) 

GO  TO  666 

27  DO  28  1=2, NCOEF 

28  COEF!I)=AN!I-1 )*GAMMA**!I-1 ) 

C 

C  COMPUTE  ROOTS  OF  POLYNOMIAL  WITH  IMSL  ROUTINE  TO 
C  OBTAIN  THE  EIGENVALUES 

C 

666  CALL  ZCPOLY(COEF, NDEG, LAMBDA, IER) 

C 

C  NOW  ORDER  THE  EIGENVALUES  3Y  SIZE 
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IF(LTEST.EQ.O)  WRITE(8,89) 

1=1 

DO  70  11=2, NDEG 

SIZE:  REAL! LAMBDA! I ) ) **2+AIMAG(LAMBDA( I ) )  **2 

K  =  I 

DO  75  J  =  1 1 ,  NDEG 

SIZE1 = REAL! LAMBDA! J ) )  **2<-AI  MAG!  LAMBDA!  J) ) **2 
IF(SIZE1 .LT.SIZE)  GO  TO  75 
K= J 

SIZE=SIZE1 
75  CONTINUE 

CDUM: LAMBDA ( I ) 

LAMBDA ( I ) = LAMBDA ( K ) 

LAMBDA  ( K )  =  CDUM. 

CL( I) rLAMBDAC I) 

.  EVPH: ATAH2! AIMAG! CL! I) ) , REAL! CL! I ) ) ) *1 80 ./PI 
SMA=REAL(CL( I) ) * *  2  +  AIMAG( CL( I ) ) **2 
SMAG=SQRT(SMA) 

IFtLTEST. EQ. 0)  V/HITEC  8 ,333)  1, LAMBDA ( I) ,SMAG, EVPH 
333  FORMAT! IX, 1 10, 4(G14.7, IX) ,/) 

1  =  11 

70  CONTINUE 

EVPH: ATAN2! AIMAG (LAMBDA (NDEG) ) , REAL( LAMBDA! NDEG) ) ) *  1 80 . /PI 
SMA= REAL! LAMBDA! NDEG) )  **24- AIM  AG!  LAMBDA!  NDEG)  )*  *2 
SMAG=SQRT(S"A) 

IF(LTEST.EQ.O)  WRITE! 8 , 333) NDEG, LAMBDA! NDEG) , SMAG, EVPH 

C 

C  NOW  CALCULATE  THE  CONSTANTS  FOR  THE  FUNCTION  SUM  FOR  A 
C  PRTICULAR  MODE.  LOOP  THEN  TO  CALCULATE  THE  FIELD  AI  A  SELECTED 
C  NUMBER  OF  POINTS  FROM  ZERO  TO  ONE 

C 

45  X=0. 

BRIGHTrO. 

WRITE(8 ,997 ) 

997  FORMAT! IX, ‘INPUT  1  TO  CALC  FIELDS,  0  TO  DO  NEW  CAVITY:*,/) 

READ  *  , MTEST2 
WRITE! 8 ,979 ) MTEST2 
IFCMTEST2 . EG . 0)  GO  TO  777 
WRITE! 8 ,995 ) 

995  FORMAT! IX,* INPUT  DESIRED  MODE  NUMBER:*,/) 

READ  * , MODE 
WRITE! 8, 979) MODE 
LABEL!  15)  =  1 0H  MODE  II 
ENCODE! 10, 937 ,LA3EL! 16) )MODE 
987  FORMAT! 12, 8X) 

NPOINT=0 

NLsNBIG 

ROOT=L AKBDA!  MODE ) 

992  F0RMAT(1X,I5,1X,*P0I»TS  WILL  BE  EVALUATED  FOR  MODE  *,2G14.7,/) 

C 

C  CONSIDER  I/O  OPTIONS  AND  CALCULATE  THE  CONSTANTS  FOR  EITHER 
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PARITY  CHOICE. 

IFCKTEST1 .EQ. 1 )  CO  TO  40 
DO  30  1=1 , ML 
RINDEXC I ) = I 

CONSTC I ) =H* ( RCOT-GAMMA) /GAMMA/ AW ( NDEG) * (  ROOT/GAMMA) ** ( NL-I) 

CO  TO  29 
DO  41  1=1,  KL 
RINDEXC I ) = I 

CONSTC I ) = (GAMMA/ ROOT )  ** I 
VRITEC8.982) 

?  FORMATC IX, *  INPUT  0  TO  CALC  INTENSITIES  OVER  EXPANDED  RANGE:*,/) 
READ  * , JTEST 
WRITEC  8 ,979 ) J TEST 
ENCODEC  1 0 , 994 , LABELC  2 ) ) NEC 
ENCODEC 10,994, LABELC 4 ) ) NAG 
ENCODEC  10,994,  LA  EEL,  (7  )  )  REALC  SCOT) 

ENCODEC  10,994,  LABELC  C  ) )  AIMAGC  ROOT) 

IF( JTEST. NE.O)  GO  TO  46 
LAB  EL  ( 9 )  =  1 0HM  IS  ?.  OR 
LAB EL ( 1 0 )  =  1 0H PLANE 
LABELC 11)  =  lOiiSCALED 
IF(IGAINQ. EQ. 1 )  LA3ELC 1 1 ) = 10H 
LABELC  12)  =  1Ci!  INTENSITY 
LABEL (  1 7 )  "  1 0 HAP ? SOX  if  2 

CALL  ALLINTC  MAG ,  MS  JBN  ,  MS'JPN ,  CONST ,  T,  N3IG ,  MTEST1  ,  ROOT , LABEL ,  H , 

1  GAMMA, IGAINQ) 

WRITEC  8,981) 

I  FORMATC IX, *  INPUT  0  TO  CONTINUE  WITH  OTHER  I/O  OPTIONS:*,/) 

READ  * , JTEST1 

WRITEC  8 ,979 ) JTEST1 

IFC JTEST 1 .NE.O)  GO  TO  45 

LABELC  1 7  )  =  1 CKAPPROX  It  1 

WRITEC8.980) 

)  FORMATC IX,* INPUT  ¥  PNTS  FROM  0-1  AND  0-1  TO  PLOT  OR  PRINT:*,/) 
READ  * , INCX , MTEST3 
WRITEC 8 ,978 ) INCX , MTEST3 
WRITEC8.992) 

CALCULATE  THE  FIELD  AT  VARIOUS  X  VALUES  USING  THE 
CONSTANTS  JUST  CALCULATED: 

FOR  EVEN,  F( X) =  1  *  3UM( H ( X ) ) 

FOR  ODD,  F( X ) =  SUM(HCX)),  WHERE  THE  H(X)'S  ARE  THE  ONES 
DERIVED  FOR  EACH  CASE 
N  POI  NT=  N  PC  I T+ 1 
STOREXC  NPOINT) =X 
SIG=CMPLX ( 0 . ,0.) 

DO  32  1=1 ,NL 

BN1=RTEYE52*SQRTCPI*T/MSU3N(I)) 

BN2  =  -T*  EYE/ MSUBNC I ) 

BN3=1 .-X/MSUPNC 1+1 ) 


BN4  =  1  .+X/MSUPNC  1  +  1 ) 

HNX=CCEXPCBr!2*5NR**2)/Bfn+CEXP(nN2*BN4**2)/B!!4)/BN1 

funvalci,:jpoint)  =  eealc-h;;x)**2+;imagC-h.hx)**2 

IFC  MTEST 1 . EQ . 0 )  GO  TO  32 

HNX=CCEXP(3N2*BN3**2)/B::3-CEXPCBM2*BN4*»2)/BN4)/3N1 
FUNVALC I , NPCI NT ) =  P.EALC  -HNX)  *  *2+AIMAG(  -HNX)**2 

32  sig=sig+co;;st(  i 

IFCMTEST1 . EQ. 1 )  GO  TO  33 
FIELDXC  N POINT ) =H+SIG 
GO  TO  34 

33  FIELDX(JIPOINT) -SIG 

34  INTENSC NPOINT) =REALC FIELDXC NPOINT) ) **2+AIMAG(FIELDX( NPOINT) ) *«2 
IFC INTENSC  NPOINT) . GT. BRIGHT )  5EIGHT=INT£NS( NPOINT) 

IFCMTEST3 . EQ. 0 )  GO  TO  <5 
WRITE ( 8 , 87) X, FIELDXC  NPCI NT) 

WRITE(  8,86)  INTENSC NPOINT) 

86  formatcix, ‘intensity  =  *, 014.7,/) 

87  FORMATC 1  X ,  ‘X  -  *  ,  G 1 4 . 7  ,  *  FIELD  =  ‘,2014.7) 

35  x=x+i./i;icx 

IF(  NPOINT.  LT.  INCX)  GO  TO  31 
IF(MTEST3 . EQ . 1 )  GO  TO  777 
WRITEC  8 ,991) 

991  FORMATC IX , *  TYPE  ZERO  TO  PLOT  CONSTANTS  VS  N:*,/) 

READ  * , MTEST C 
WRITEC  8, 979  )M.7ESTC 
IFCMTESTC. NE.O)  GO  TO  38 
LABELC 9 )  =  1 0HCONSTA.'iT  l) 

LABELC  10)  =  1C!i 
LABELC  IDslOH  MOD (  CONS 
LABELC 1 2 ) s 1 CHT  ANT) **2 
DO  42  1=1,  fIL 

42  PLOCOHC I )  sBEALC  CONSTC I ) )  “2+AIMAC  ( COfiSTC  I) )  “2 
CALL  KGRAPHCRINDEX.PLCCON.NBIG.LAEEL,  1  ,-1  ,  1  1) 

WRITEC 8 ,984) MODE 

984  F0RMATC1X, ‘COMPLETED  PLOT  OF  CONSTANTS,  MODE  =*,I2,/) 

38  WRITE ( 8 , 990 ) 

990  FORMATC 1 X , *TY  PE  INDEX  CF  FUNCTION  TO  PLOT  OR  0  TO  CONTINUE:*,/) 
READ  *, INDEX 
WRITE (8, 979) INDEX 
LABELC 9 )=10H MIRROR 
LABELC 10) = 10 H PLANE 
IFC INDEX . EO . 0 )  GO  TO  36 
DO  43  1=1 , INCX 

43  PLOFUNC I) sFUNVALC INDEX , I) 

LABELC  1  1  )=10HMOD(FUN,IN 
ENCODF.C  1 0 ,989  , LABELC  1 2 ) )  INDEX 

989  FORMATC "DEX=" , 12, " ) *  *  2 " ) 

CALL  HG R AF H C  STO R E X , P LO FU N , I NCX , L ABEL , 1 , 0 , 0 ) 

WRITEC8,986)INDEX 

986  FORMATC IX, ‘COMPLETED  PLOT  OF  FUNCTION,  INDEX  =*,I2,/) 

GO  TO  38 
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36  WRITE(8,988) 

988  FORMAT( 1 X , ‘TYPE  ZERO  TO  PLOT  INTENSITY:*,/) 

READ  * , ICONT 
WRITE( 8 ,979 ) ICONT 
IF( ICONT. NE.O)  GO  TO  45 
LABEL ( 1 1 ) : 1 OHSCALED 
IFdGAIHQ.EQ.  1 )  LABEL! 11 ): 1 0H 
LABEL! 1 2 ) : 1 0H  INTENSITY 
IF( IGAINQ. EO . 1 )  BRICHTsl. 

DO  37  1= 1 , INCX 

37  INTENS(I)  =  INTENS(  D/BRIGHT 

CALL  HGRAPHC  STOREX,  INTEL'S,  INCX,  LABEL,  1 ,0,0) 

WRITE( 8 ,985 ) 

GO  TO  45 

978  FORMAT  (IX,  *  IN' PUT  VALUES  ARE  :  *,215,/) 

979  FORMAT! IX, *  INPUT  VALUE  IS  :  *,I5,/) 

985  FORMAT( IX, COMPLETED  PLOT  OF  NORMALIZED  INTENSITY.*,/) 

998  FORM AT( 1 X , *  REVISE  PARAMETERS  SO  NDEG  50*) 

88  FORMAT( 1 CX , *  MAG  =  * , F6 . 2 , 5X , *N£Q  :  *,F6.2,/) 

89  FORMAT! 9X , *  I  * , 2X , *LAM3DA( REAL) * , 2X , ‘LAMBDA ( IMAG) * , 6X , 

1  *EVMAG*,11X,*EVPH*,/) 

888  CALL  EXIT 
END 

SUBROUTINE  ALLINTC MAG, MSU3N , MSUPH , CONST ,T , NBIG , MTEST1 , ROOT 
1  , LABEL, H, GAMMA, IGAINQ) 

DIMENSION  LABEL! 17), XSAVE! 2000) 

REAL  MS U B t.' ( 5 1 )  ,  MSUPN(  5  1 )  ,  INARG 1  ,  JUARG2 ,  INARG 3  , INARG4  ,  INARG5 
REAL  INARG6  ,  INARG7  ,  INARGd ,  MAG ,  If.'TENSC  2000 ) 

REAL  MINV 

COMPLEX  APART1 , APART2 , B PARTI , BPART2 , ALLFUN , CONST! 5 1 ), ROOT , EYE 
COMPLEX  AFUN ,  BFUN ,  SP.'.'TC ,  SPNTD  ,  EVENX ,  OUTCON ,  FRESL,  CONSTA ,  CONSTB 
COMPLEX  EYEFAC , SPCON 

C 

c 

C  THIS  SUBROUTINE  FOLLOWS  PROGRAM  BARC  AND  COMPUTES  BEAM  INTENSITITE 
C  IN  THE  OUTPUT  PLANE  FPOM  THE  CPTIC  AXIS  TO  SOME  DESIRED  OUTER  POIN 
C  OUTER  POINT  AND  i  INTERMEDIATE  POINTS  FOR  EVALUATION  ARE  INPUT 

C  WHILE  ALL  OTHER  REQUIRED  QUANTITIES  ARE  CARRIED  THROUGH  III  THE 

C  ARGUMENT  LIST  AS  FOLLOWS: 

C 

C  MAG:  CAVITY  MAGNIFICATION 

C  MSUBM:  ARRAY  FOR  PARTIAL  SUMS  OF  INVERSE  POWERS  OF  MAG 
C  HSUPN:  ARRAY  FOR  MAG  TO  SOME  POKER 

C  CONST=  ARRAY  CF  CONSTANTS  IN  THE  ASYMPTOTIC  SERIES 

C  T:  QUANTITY  DEFINED  IN  BAP.C  PER  HORWITZ 
C  NBIG:  0  TERMS  IN  THE  SERIES 
C  MTEST1 :  PARITY  DESIGNATOR 
C  ROOT:  MODE  EIGENVALUE 
C  LABEL:  PLOT  LABELING  ARRAY 
C 
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BRICHTrO. 

PI=2.«ASINC1 .) 

EYE=CMPLX( 0 . ,  1  . ) 

EYEFAC=C1 .-EYE)/2. 

WRITE( 8 , 900 ) 

900  FORMAT! INI , ‘ENTERING  EXTENDED  RANGE  INTENSITY  SUBROUTINE.*,/) 

DO  10  1=1,51 

10  MSUPfK  I)  =  f!AG*MSUPNC  I) 

NPOINT= 1 
X=0. 

WRITEC  8 ,901) 

901  FORMAT( 1 X , *  INPUT  MIN  AND  MAX  X  VALUES  AND  #  POINTS  BETWEEN:  *,/) 
READ  * , XMIH , XMAX , INCX 

X=XMIN 

50  ALLFUN=(0. ,0. ) 

X0MAG=X/MAG 
DO  310  1=1 ,NBIG 
MINVr 1 ./MSUPN( I ) 

EVENX=(0. ,0. ) 

SPNTD=(0. ,0. ) 

SPNTC= ( 0 . ,0. ) 

CONST A=- CONST! I) *GAMMA 
CONSTBrCONSTA 

IF(MTEST1 .EQ. 1)  CONSTBr-CONSTA 
P2PRYMr  2 .  *  C  1  .  +  1  .  /  MSU  PNC  2*  I )  /  MS'JBNC  I ) ) 

STAPHAz ( MINV/MSU3NC  I )+ XCMAG )/(  ,5SP2?RYM) 

STAPHBz  C -MI N7/MSU3N ( I ) +XOMAG) / ( . 5*  P2PRYK) 

INARG1  =  C 1 . -XCMAG) *  *?  +  ( (  1  .  -  MI  NV  )  *  *2/ MSU BN  CD) 

INARG2= ( ( 1 . -XOMAG-C  MINV-MINV* *2 )/MSUSNC I) ) **2 ) /( . 5  * P2PRYM) 

AARG1 = I NARG2-IN ARG 1 

APART1  =  CEXP(  EYt*T* A AF.G  1 ) / ( 1  . -MINV )*(-CO!;STA) 

INARG3  =  (  -  1  .  -  XCMAG )  *  *2*C  1  .  +MINV  )  *  *  2/MSUB(.'C  I ) 

INARG4=C(-1 . -XOMAG-C MINV  +  MINV  **2 )/ ( MSU3N( I ) ) ) **2)/( .5*P2PRYM) 
AARG2  =  If'ARC-il  - 1 NARG3 

A  P  A  R  T  2  =  C  E  X  P  (  E  Y  E  *  T  *  A  A  R  0  2 )  *  C  -  C  0  N  S  T  A )  /  ( 1  .  +  M I  f  I V ) 

INARG5  =  (  1  .-XCMAG)  **2  +  C  1  .  +  MIN7  )  *  * 2/MSUBNC I ) 

INARG6  =  ( 1  .  -  XCMA S  +  C  N. I NV  +  "  CNV  *  *2  )  /  MSUBNC I )  )  **2/  (  .  5*  P2PRYM) 

BARG 1  =  1 NARG6 - 1 N A  EG  5 

BPARTIzCEXPCSAP.OI *  EYF.*T>  * C -CONST? )/ ( 1 ,+MINV) 

INARG7=  ( 1  . +XOMAG)  **?+(’  .  -  MINV )  *4-2/MSU5N(  I ) 

INARG8  =  ( -1 .  -  XOMAG+C  UIN7-.MINV**2 ) /  MS'JENC  I) )  **2/  ( . 5* P2PRYM) 
BARG2=INARG8-INARG7 

BPART2  =  CEXP(BAr<G2*EYE*  r)  *(-C0N3TB)/(  1  .-MINV) 

OUTCON  =  SQRT( MSU3N ( I ) / PI/T/P2PR Y!‘. ) /  2 . / ROOT 
SPCONzSQP.T (  MS'JBNC  I )  *2 .  /  P2PRYM/ PI/T) /  2 . / FOOT 

FRSPOS=SQRT C  T/P 1/ P2PRYM) * ( 2. * ( 1 . -XOMAG ) -2 . * C 1 . -Ill NV ) /KSUPNC I) / 

1  KSUBKCD) 

FRSNEG=SQRT(T/PI/P2PRYK)*C2.*C-1.-X0MAG)-2.*(1 . * HI NV )/MSUPN( I ) / 

1  MSUBNCI)) 

IFCSTAPHA.GE. -1 . . AND. STAPMA.LE. 1 . )  GO  TO  130 

AFUN=APART1*CFRESL(FRSP0S)-EYEFAC)-APART2*(FRESL(F RSNEG ) -EYEFAC ) 
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IFCSTAPHA.LT. 1 . )  GO  TO  200 

AFUNsAPARTI *  C  FRESLC FRSPOS) +EYEFAC )-APART2* ( FRESLC FRSNEG )+EYEFAC) 
GO  TO  200 

130  AFUNsAPARTI * C  FRESLC FRSPOS ) -EYEFAC )- APART2* ( FRESLC  FRSNEG 5 +EYEFAC) 
SPNTCrSPCON*(-CONSTA)/(  1  . -STAPHA/ KSUPNC I )  ) *CEXPC -EYE* C PI/4 . +T* 

1  (  ( STAPHA-XCMAG ) **2+( 1  . -STAPHA/MSUPNC  I)  )  **2/KSUSNC  I )  )  )  ) 

200  FRSPOSsSQP.T  CT/PI/P2PRYM)*C2.*C1  .  -XOKAG)  +2.  *(  1 . +MINV)/MSUPHC I)/ 

1  MSUBN(I)) 

FRSNEG=SQRTCT/PI/P2PRYM)*C2.«C-1 ,-XOMAG) +2. * ( 1 .-MINV )/MSUPNC I )/ 

1  MSUBNCI)) 

if(staphb.ge.-i..and.staphb.le.i.)  GO  TO  140 

BFUNsB PARTI  * ( FRESL C FRSPOS) -EYEFAC ) -3PART2* ( FRESLC  FRSNEG ) -EYEFAC) 
IFCSTAPHB.LT. 1 . )  GO  TO  300 

BFUNsB PARTI *( FRESL ( FRSPOS) +EYEFAC )-BPART2* ( FRESLC FRSNEG ) +EYEFAC ) 
GO  TO  300 

140  BFUN=BPART1*C FRESLC FRSPOS) -EYEFAC )-BPART2*C FRESLC FRSNEG )+EYEFAC) 

SPNTDsSPCON* C -CONSTB)/(  1 .  +STAPHB/MSUPNC  I ) )  *CEXPC -EYE*  (  PI/ 4 .  +T* 

1  ( ( STAPHB-XOXAG ) *  *2+ ( 1 . +STAPH3/MSUPNC I ) ) *  *2/KSUBN ( I ) )  )  ) 

300  ALLFUNsOUTCON*(  AFUN+EFUN)+SPNTC+SPNTD+ALLFUN 

310  CONTINUE 

802  FORMATC IX ,2G14 .7) 

EARG1sSQRT(T/2./PI) *2.*(1 .-XOMAG) 

EARG2  =  SQP,T(  T/2 .  /  PI)  *2  .  *  ( - 1  ,-XCMAG) 

EVENXsCSQRTC EYE/2 . }/ ROOT* ( FRESLC  EARG1 ) -FRESLC EARG2) ) 

IFC X/MAG. GE.  - 1  . .  AND .  X/MAG.  LE.  1 . )  EVENXs EVE.NX-CSQRTC  EYE/2 . ) /ROOT* 
1  (1.-EYE)+CEXP(-EYE*PI/4.)/R00T*CSQRT(EYE) 

EVENXs EVENX*H*GAMM A 

IFCMTEST1 .EQ.O)  AU.F UH= ALLFUN+EVENX 

WRITEC  8 , 802 ) ALLFUN 

INTENSC  NPOINT!  sAIMAGC  ALLFUN) * *2  + REALC  ALLFUN)  **2 
xsavecnpoi:;t)sx 

WRITEC 8, 800) INTENSC MPOINT) ,XSAVE( NPOINT) 

800  FORMATC 1 X , 2G 1 4 . 7 ) 

IFC INTENSC  NPO I  NT) . GT . E RIGHT)  BRIGHTsINTENSC NPOINT) 

IFCIGAINQ. EQ. 1 )  BRIGHTsl. 

IF(X.GT.XMAX)  GO  TO  500 
XsX+1 ,/INCX 
NPOINTsNPOINT+1 
GO  TO  50 

500  DO  510  1=1  .NPOINT 

510  INTENSC I) = INTENSC I) /BRIGHT 

CALL  HGRAPHC  XS A VE ,  INTEL'S ,  NPOINT  .LABEL ,  1  ,0,0) 

DO  600  1:1 ,51 

600  MSUPNC I )  sMSUPNC I) /MAG 
WRITEC  8,904) 

904  FORMATC IX, ‘COMPLETED  CALCULATION  AND  PLOT,  EXTENDED.*,/) 

RETURN 

END 

SUBROUTINE  HGR APHC X , Y , N, ID , NO , NP , NS) 

DIMENSION  X( 1 )  ,  YC 1  )  ,  ID( 1 )  $  IF  (N0.EQ.2)  CALL  PLOTC -1 . 85 , 2 . 1 0 , -3 
IF  (N0.EQ.2)  GO  TO  30  $  IF  CNO.LT.O)  GO  TO  10 
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CALL  SCALECX.7. ,N,1)  $  CALL  SC ALEC Y , 5 . , N , 1 ) 

10  CALL  PL0T(0. , 11 .  ,2)  $  CALL  PL0T( 8. 5 , 1 1 . ,2) 

CALL  PL0T( 8 . 5 i 0 .  ,2)  S  CALL  PL0T( 0. , 0 . , 2 ) 

CALL  PLOTC 1.35,1.35,-3)  $  CALL  PL0T(0. ,8. 30, -2) 

IF( ID( 1 ) . EQ.999 )  CO  TO  25 

CALL  PLOT ( . 1 , - . 1 , -3 )  $  CALL  PLOTC 0 . , -2 . , -2) 

DO  20  1=1  ,7,2 

20  CALL  SYMBOLC  ( 1+ 1 . 5) * . 1 , . 3 , . 07 , ID( I ) , 90 . , 20) 

CALL  PLOTCO. ,0. ,3)  $  CALL  PLOTC 1 . , 0 . ,2 ) 

CALL  PLOTC 1 . ,2. ,2)  $  CALL  PLOTC 0 . ,2 . , -2) 

CALL  PLOTC-. 1 , . 1 ,-3) 

25  CALL  PLOTC  5.8,0. ,-2) 

CALL  PLOTCO. ,-8.30, -2)  $  CALL  PLOTC -5 . 8 , 0 . , -2) 

CALL  SYMBOLC .5,-.2, .1  ,IDC13) ,0. ,50)  $  CALL  PLOTC 5 . 3 , . 75 , -3 ) 

CALL  AXISCO. ,0. , IDC 9) ,-20 ,7 . ,90. ,X< N+1 ) ,X(ti+2) ) 

CALL  AXISCO. , 0.  ,  IDC  1  1 )  ,20,5.  ,180.  ,  Y( .'!+ 1 ) ,  Y(  M+2) ) 

30  Y C  N+2)  =- Y(  !l+2 )  t  CALL  LINEC  Y ,  X ,  N ,  1  ,  TIP ,  NS) 

Y( N+2 ) =- Y( N+2 )  3  CALL  PLOTC 1 . 85 , -2 . 1 0 , -3) 

RETURN  S  END 

SUBROUTINE  AX  ISC  XO  ,  YO , L, NC , RL , ANG , RMIN , DR) 

DIMENSION  LCD  $  A=  ANG  *3  .  1 « 1  59/ 1 80.  $  DX= . 1*C0S( A)  $  DY=.1*SIJJ(A 
ICrISIGHCl ,NC)  $  MNCsIABSCHC)  S  R:.1  $  Nsl  $  X=XO  $  Y=Y0$ 

10  CALL  PLOTC  X , Y , 3'  3  XsX+DX  3  Y  =  Y+DY  $  CALL  PL0T(X,Y,2) 

CALL  PLOTC X-. 21 »DY*IC, Y+.21 «DX*IC,2) 

IFCN.EQ.5)  CALL  PLOTC  X- .  *12*D :  *IC ,  i  +  .  U2*DX*IC  ,2) 

IF(N.EQ.IO)  CALL  PLOTC X- . 70*DY*IC , Y+ . 70 *DX*IC , 2) 

NrMODCN, 10)+1  $  R=R+.1  3  IFCR.LT.RL)  GO  TO  10 
A=ANG-CIC+1 )*45.  3  DX= 1 0 . *DX  $  DY=10.*DY 
C=-. 175+.  1 25  * IC  $  D= . 19  + . 35* IC 

-  X=XO+C*DX-D*DY  $  Y=YO+C*DY+D*DX 

RsAMAXI  (  ABSC  FMIfl)  ,  ASSC  3t'I!.'+DR*RL) )  $  RsALCGIOCR) 

IRsINTCABSC  R) )  3  IFCR.LT.O.)  IH=-CIR+1)  $  IR  =  IR-M0D(IR,3> 

R1 sRMIN/1 0 . ** IR  $  DR  1 =DR/ 1 0 . *  * IR  $  R  =  0. 

20  ENC0DEC7 , 101 , S) R 1  3  CALL  SYMBOLC X , Y , . 07 , S, A, 7 )  $  R1=R1+DR1 
X=X+DX  $  YsY+DY  3  R=R+1.  3  IFCR.LE.RL)  GO  TO  20 
R=  (RL- .  1  *  NNC)  /2 .  $  Cr.U.5*IC 

X=X0+R*DX-C*DY  $  Y=YO+R*DY+C*DX 

CALL  SYMBOLC X , Y , . 1  , L , ANG , NNC)  3  IFCIR.EQ.O)  RETURN 
ENCODEC 5 , 1 02 , S)  $  CALL  SYMBOLC 999 . ,999 .,. 1 0 , S, ANG , 5) 

CALL  WHEREC X , Y , A) 

ENCODEC 3 , 103 , S)  IR  $  CALL  SYMBOLCX.Y, .07,S,ANG,3) 

101  FORMATC  F7 . 2 ) 

102  FORMATC  5H  *10) 

103  FORMATC 13) 

RETURN  $  END 

C  s  SUBROUTINE  SCALEC DATA , LENGTH , N , K)  r 

C  =  = 

C  =  REAL  DATA  =  N+2  DIMENSIONED  ARRAY  OF  DATA  TO  BE  SCALED  = 

C  =  INTEGER  N  =  NUMBER  OF  DATA  POINTS  = 

C  =  REAL  LENGTH  =  LENGTH  OF  THE  PLOT  AXIS  (E.G.  IN  INCHES)  = 

C  =  INTEGER  K  =  UNUSED  PARAMETER  INCLUDED  FOR  COMPATIBILITY  = 
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s  WITH  THE  EQUIVALENT  CALCOHP  SUBROUTINE  s 
=  s 
=  THE  FOLLOWING  VALUES  ARE  RETURNED:  s 

s  DATACN+1)  =  ADJUSTED  DATA  MINIMUM  = 
=  DATA( N+2)  =  "NICE"  SCALE  FACTOR  IN  DATA  UNITS  = 
=  PER  LENGTH  UNIT  (E.G.  VOLTS/INCH)  = 


SUBROUTINE  SCALE( DATA , LENGTH, N , K) 

REAL  DATA(N) ,  LENGTH,  SF(5) 

DATA  SF  /I . ,  2. ,  2.5,  5.  ,  10.  / 

COMPUTE  THE  RAW  SCALE  FACTOR 

DMIN=DMAX  =  DATA(1  ) 

DO  10  1=1, N 

IF(DATA(I)  .LT.  DMIN)  DM IN  =  DATA(I) 
IF  (DATA(I).  GT.  DMAX)  DMAX  =  DATA(I) 
0  CONTINUE 


EXCLUDE  TRIVIAL  ERROR  CASES 
DATA( H+ 1 )  =  DMIN 
DATAOI+2)  =  1.0 

IF  (LENGTH  .LE.  0.0  .OR.  DMAX  .EQ.  DMIN  )  RETURN 


RAWSF  =  (DMAX  -  DMIN)  /  LENGTH 

C  RAWSF  =  SFMANT  *  10.  **  SFEXP,  WHERE  1  .LE.  SFMANT  .LT.  10 

SFEXP  =  AINTC  AL0G10C  RAWSF  )  ) 

IF  (  RAWSF  .LT.  1.0  )  SFEXP  =  SFEXP  -  1.0 
SFMANT  =  RAWSF  *  10,0  **  (-SFEXP) 

C  LOCATE  NEXT  LARGER  "NICE"  SCALE  FACTOR 

DO  20  1=1,5 

20  IF  (  SF( I )  .GT.  SFMANT  )  GO  TO  30 

PRINT*,"  SCALE:  SCALE  FACTOR  ERROR  ..."  $  RETURN 

30  SFNICE  =  SF( I )  *  10.0  «*  SFEXP 

C  COMPUTE  ADJUSTED  DATA  MINIMUM 

AD J MIN  =  AINT  (  DMIN  /  SFNICE  )  *  SFNICE 

IF  (  ADJMIN  .GT.  DMIN  )  ADJMIN  =  ADJMIN  -  SFNICE 

IF  (  (DMAX  -  ADJMIN)  /  SFNICE  . LE .  LENGTH  )  GO  TO  HO 
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C  NEED  TO  USE  THE  NEXT  LARGER  SCALE  FACTOR 

IF  (  I  .LT.  5  )  SFNICE  =  SPCI+1)  *  10.0  **  SFEXP 
IF  (  I  .EQ.  5  )  SFNICE  =  20.0  *  10.0  **  SFEXP 

ADJMIN  =  AINT  (  DMIN  /  SFNICE  )  *  SFNICE 
IF  (  ADJMIN  .GT.  DMIN)  ADJMIN  =  ADJMIN  -  SFNICE 

40  CONTINUE 

DATACN+1)  =  ADJMIN 

DATA( N+2)  r  SFNICE 

RETURN 

END 

COMPLEX  FUNCTION  CERF(ZZ) 

COMPLEX  ZZ,Z,A,A1 ,A2,B,B1 ,S2,F,F1 
Z=ZZ 

IF(CABS(Z) . CE. 3 . 5G0T030 
J=0. 

A=Z 

BsZ 

10  J=J+1 

B=-Z*Z«CMPLX(FL0AT(2*J-1) ,0. )«B 

B=B/CMPLX( FLOATC  J  ) , 0. ) /CKPLXC  FLOATC 2* J  +  1 ) , 0. ) 

AsA+B 

IF(J.GE.1000)GOT050 

IF(CABS(B/ A) . GE . ( 1 . E- 1 0 ) )  GO  TO  10 

CERF=(1. 128379167,0. )*A 

RETURN 

30  IF(REAL(ZZ).LT.O.)Zsr-ZZ 
A2=C1.,0.) 

B2=Z 

F2=A2/B2 
A1  sZ 

BUZ»Z  +  (0.5,0.) 

F1=A1/B1 

J=1 

40  J=J+1 

A=Z*A1+CMPLX( FLOATC J)/2. ,0.  }*A2 
B=Z*B1+CMPLX(FL0AT( J)/2. ,0. )*B2 

F=A/B 

IF(J.GT.1O0O)GOTO  50 

IFCCABSC ( F-F1 )/F).LT.(1.E-l0))G0T060 

A2rA1 

B2=B1 

AlsA 

B1=B 

F1  =  F 

GOTOUO 

50  WRITE(8,99) 
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99  FORMAK  "  ERROR  FUNCTION  ROUTINE  DID  NOT  CONVERGE  ") 
IERsI 
RETURN 

60  F1s(0.5,0.)*CEXP(-Z*Z)*F 
CERF=1 . 128379 1 67*F1 
CERFrl .-CERF 

IF(REAL(ZZ).LT.O.)  CERF=-CERF 
70  RETURN 
END 

COMPLEX  FUNCTION  FRESL(X) 

COMPLEX  EYE , Z , CERF 

EYE= (0 . , 1 . )  $  PI=2.»ASIN(1 .) 

Z=SQRT(PI)*X*(1 . -EYE) /2. 

FRESLr ( 1.+EYE)/2.*CERF(Z) 

FRESL=CONJG( FRESL) 

RETURN  $  END 
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